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SUMMARY 


The numerical solution of the Euler or Navier-Stokes equations by Lagrangian “vortex 
methods” is discussed. The mathematical background is presented in an elementary fash- 
ion and includes the relationship with traditional point-vortex studies, the convergence to 
smooth solutions of the Euler equations, and the essential differences between two- and 
three-dimensional cases. The difficulties in extending the method to viscous or compress- 
ible flows are explained. 

The overlap with the excellent review articles available is kept to a minimum and more 
emphasis is placed on the author’s area of expertise, namely two-dimensional flows around 
bluff bodies. When solid walls are present, complete mathematical results are not available 
and one must adopt a more heuristic attitude. The imposition of inviscid and viscous 
boundary conditions without conformal mappings or image vortices and the creation of 
vorticity along solid walls are examined in detail. Methods for boundary-layer treatment 
and the question of the Kutta condition are discussed. 

Practical aspects and tips helpful in creating a method that really works are explained. 
The topics include the robustness of the method and the assessment of accuracy, vortex- 
core profiles, time-marching schemes, numerical dissipation, and efficient programming. 
Operation counts for unbounded and periodic flows are given, and two algorithms designed 
to speed up the calculations are described. 

Calculations of flows past streamlined or bluff bodies are used as examples when appro- 
priate. These include curved mixing layers, the starting vortex and the dynamic stall of an 
airfoil, rotating stall in a two-dimensional cascade, a multi-element airfoil, and an attempt 
at predicting the drag crisis of a circular cylinder. 


1 


1. INTRODUCTION 


Vortex methods in general were thoroughly reviewed by Leonard (1980, 1985), and it 
would be difficult to improve on these articles. In the present notes we intend to cover the 
basics of vortex methods and then to discuss in more depth two subjects: the interaction 
with solid walls, and the practical aspects of programming and using vortex methods. 
These subjects were of less interest to Leonard, but are crucial in engineering applications. 
Therefore we attempt to present useful and sometimes original material in these areas. 
We shall also mention some recent contributions, published after Leonard’s articles. The 
sections which contain advanced material and could be omitted for a first reading are 
indicated by stars.* 


1.1. Example: Flow Past a Multi- Element Airfoil 

We propose to start by showing a calculation of a separated flow via the vortex method to 
illustrate its main features and, we hope, render it attractive to the reader. We shall often 
describe the method by comparison with grid-based (finite-difference or finite-element) 
methods, since these are more widespread. 

In 1986 the Boeing Commercial Airplane Company contacted the author, asking for 
computations of the flow past an airfoil in landing configuration, with a leading-edge slat 
and a double-slotted trailing-edge flap. Boeing was especially interested in the variation 
of the lift, with Reynolds number. The configuration was a severe test for numerical meth- 
ods due to the complex, multiply connected shape and the importance of viscous effects, 
including separation. 

Results were obtained in less than a week. In spite of the many concave and convex 
sharp corners present, no smoothing or other alteration of the shape was necessary. This 
illustrates the first advantage of the vortex method: it is grid-free. With a finite- element or 
especially a finite-difference method the grid generation for such a shape would be difficult 
and time consuming, and one may have to settle for rather distorted grids which degrade 
the accuracy. With the vortex method the only precaution taken was to make the spacing 
between vortices and the time step small enough to accomodate the narrow gaps between 
elements 2 and 3, and 3 and 4 of the airfoil. 

Figure 1 shows the airfoil, the vortices (more accurately, the centers of the vortex blobs), 
streamlines, and the individual (hollow arrowheads) and overall (filled arrowhead) force 
vectors at a given time. Observe the recirculating bubble behind the slat, the acceleration 
of the fluid through the slots, the separation at sharp corners and on the top surface 
of element 4, and the irregular motion in the wake. The computed forces were in good 
agreement with experiment. The figure illustrates the Lagrangian character of the method: 
the solution procedure consists in tracking vortices which move with the fluid, rather than 
updating quantities at fixed grid points. It also illustrates the second advantage of the 
method: the vortices carry all the information and are needed only in a narrow region 
near the body and in the wake. In fact the calculation used just 1300 vortices (roughly 
the equivalent of a single 36 x 36 grid). The streamlines were computed only for display 
purposes; for the purpose of advancing the equations in time, computing the forces, etc., 
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Figure 1. Flow past a multi-element airfoil. • • • vortices; f forces; — streamlines 

one only needs to compute the velocity of the vortices themselves. 

Additional advantages of the method are first that the boundaries of the computed 
domain are at infinity which removes all of the problems associated with computations in a 
truncated domain (the effects of the boundary conditions on accuracy and on stability) and 
second that there is low numerical dispersion. The method can transport flow structures 
(groups of vortices) at any velocity without deforming them or dissipating them as a grid- 
based method may do. There is no CFL number. The question of time accuracy will be 
addressed in much more detail later. 

Among the disadvantages are the fact that the flow had to be treated as incompressible, 
although fairly high local Mach numbers may be reached in the slots even at the landing 
speed. The calculation also required the storage of a 249 x 249 full matrix and about 
5 x 10 11 floating-point operations to establish the flow, starting from rest, and generate a 
time sample long enough for averaging. Both figures are much larger than what a 36 x 36 
grid calculation would require. Thus, a vortex carries much more information than a grid 
point, but it is also much more expensive to maintain. This is because all the vortices 
interact, so that the operation count is of order N 2 , where N is the number of vortices. 
The competing methods often have operation counts of order N or iVlog(iV). 

To answer Boeing’s question, a Reynolds number dependence was indeed predicted. As 
the chord Reynolds number was changed from 2.3 x 10 6 to 1.2 x 10 7 with all the other 
parameters (physical and numerical) exactly the same, the lift coefficient at 8° incidence 
increased by about 5%. This was caused by a difference in the boundary-layer behavior 
on the upper surface of element 4. At the lower Reynolds number, the boundary layer 
always separated, whereas at the higher Reynolds number it could transition and remain 
attached part of the time, thus increasing the circulation around element 4 and therefore 
the overall lift. This Reynolds-number trend was consistent with the common wisdom. 
Ironically, Boeing subsequently informed the author that experiments had predicted an 
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opposite trend for this configuration, so that the situation was rather confused. This 
study demonstrates the power of a well-designed vortex method in dealing with complex 
geometries, and how challenging and unpredictable Reynolds-number effects can be to 
numericians, experimentalists, and designers. 


1.2. Governing Equations 

Complete discussions of these equations can be found in textbooks, e.g., Batchelor 
(1967). The fundamental equations are the incompressible Navier-Stokes equations, 


V.U = 0 

U t + U.VU = -Vp+iA7 2 U. 


(la) 

(16) 


Here U is the velocity vector, V is the gradient operator, t is the time, p is the kinematic 
pressure, and u is the kinematic viscosity. Bold-face letters like U denote vector quantities, 
and the symbol = indicates a definition. Equation (la) is the continuity equation and (lb) 
the momentum equation. One obtains the Euler equations by setting u = 0. The vorticity 
is the curl of the velocity, 

u> = V x U. (2) 

In an unbounded domain with the fluid at rest in the far field, one has the “Biot-Savart 
law”, which is the inverse of (2) by a Green’s-function approach (when (la) holds) and 
explicitly gives the velocity in terms of the vorticity: 


U(x) 


1 f u>(x') x (x - x') 

47 r J jx — x'| 8 


dx'. 


(3) 


This equation only expresses the kinematics of the flow. The integral, taken over the whole 
two- or three-dimensional (2-D, 3-D) domain, is responsible for the nonlocal interactions 
between vortices. By taking the curl of (lb) one gets the vorticity equation, 


u> t + U.Vu> = u.VU + vV 2 u>. 


(4) 


Notice that the pressure term is eliminated, but that a new term, u>.VU, appears which is 
responsible for vorticity stretching and rotation. In this term, VU can be replaced by its 
symmetric part (VU + VU T )/2 (the deformation tensor), which may be advantageous in 
some numerical methods (Cantaloube and Huberson, 1984). 

The circulation of a velocity field U along a closed contour is defined by the line integral 

T = j\l.d\. (5) 

This is an extremely useful concept because Kelvin’s theorem (Batchelor, 1967, p. 273) 
shows that in the absence of viscosity, the circulation around a contour that follows the 
fluid is conserved: DT / Dt = 0. 
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In 2-D the equations are much simpler. If the flow is in the (x,y) plane the vorticity 
only has a component in the 2 direction, u> 2 = dv/dx - du/dy (which we’ll denote simply 
by u>). The Biot-Savart law becomes 


u,x) = i e ‘ x I **■' (6) 

with e. the unit vector in the 2 direction. The vorticity equation (4) becomes 

u> t -|- U.Vuj u }. (7) 

Note that the rotation/stretching term has disappeared. One can define a scalar stream 
function V’ by U = e z xVif since U is divergence- free. Then V’ satisfies a Poisson equation, 
V 2 ^ = u >, and the Biot-Savart law is written 

v>(x) = ^ J w(x')log(|x - x'l) dx'. (8) 


1.3. Invariants of the Motion 

One important source of confidence in a numerical method is its conservation of at least 
some of the appropriate integral quantities. In the present context, the most useful form 
of these invariants is in terms of the vorticity, rather than the velocity. See the discussions 
by Batchelor (1967, p. 517), Wu and Sankar (1980) or Ting (1983); the arguments rely 
on integrations by parts but are subtle, because the velocity decays slowly (algebraically) 
at large distances even when the vorticity decays fast (exponentially) which is the usual 
situation. In 3-D the important integrals are 

(9) 

( 10 ) 

( 11 ) 

(12) 

They correspond to total vorticity, linear momentum, angular momentum, and kinetic 
energy, respectively. We assume that the vorticity decays fast enough in the far field for 
these integrals to be defined. Note that the total vorticity (9) is trivial unless u; decays 
very slowly, because if it decays any faster than r -3 , one can easily show that this integral 
is identically zero (use V.u> = 0). In 2-D the corresponding quantities are 

J w dx, (13) 


J u> dx, 

J x x u; dx, 

J x x (x x u>) dx, 

( dx dx'. 

J |x-x'| 
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(14) 


-e 2 x J u x dx, 



J u>(x).u>(x') log(|x - x'l) dx dx'. (16) 


Here, the total vorticity (13) has no reason to be zero. 

In unbounded inviscid flows, all these quantities are conserved. In unbounded viscous 
flows, only total vorticity and linear momentum (9, 10, 13, 14) are conserved. Formulas 
giving the rate of viscous decay of angular momentum and kinetic energy can be found 
(e.g., Ting, 1983). When there is a solid body in the fluid it can alter any of the integrals 
(9-16). However, the variation of (9) and (13) (with the integrals restricted to the fluid 
region) is known once the motion of the body is (Wu and Sankar, 1980). 

The invariants are often invoked when constructing or evaluating methods, but the 
actual importance of exactly conserving quantities such as (9-16) in a numerical method 
is a matter of debate. The invariants are integrals over the whole domain, and contain no 
information about the details of the solution. Many successful methods violate some of 
the invariants, especially once time-integration errors are included, and one can construct 
methods that conserve the invariants but are not very good or even are inconsistent. The 
common argument that energy conservation prevents the calculation from “blowing up” is 
rather technical, and in any case does not directly apply to vortex methods. 

The following intuitive argument seems more relevant. The conservation of an invariant 
means that over the whole domain the local errors exactly cancel each other in some 
measure (e.g., |x 2 |dx for (15)). From a statistical point of view, it is then reasonable to 
expect that over an intermediate domain, large enough to contain many vortices or grid 
point s, the local errors tend to cancel each other instead of accumulating. One then expects 
much better accuracy on intermediate scales than on small scales. This is probably true 
in many vortex calculations, especially complex ones such as in figure 1. One can tolerate 
significant errors in the irregular small-scale motion, as long as they do not affect the 
quantities of interest such as the forces and moments on the bodies. 


1.^. Constraints on Vector Fields * 

Vortex methods and related methods focus on derivatives (or combinations of deriva- 
tives) of the primary variables. Using a collection of computational elements (e.g., vortices) 
one can represent an arbitrary vector field. This can become a liability, because many of 
the vector fields involved are not arbitrary; they satisfy constraints. For instance, if the 
vector A is the gradient of a scalar, say A = V<^>, then V x A = 0. Now if one represents 
the components of A separately with two collections of computational elements there will 
in general be an inconsistency because dA x /dy may not be equal to dA y /dx , and if they do 
differ there is no scalar <f> such that d<j>/dx — A x and d<t>/dy = A y . For instance Anderson 
(1985a) proposed a Lagrangian method for the transport of a gradient and avoided this 
question by solving a Poisson equation for <f> : V 2 ^> = dA x /dx + dA y /dy , instead of the 
original equations. For the vorticity, the constraint is that it should be divergence-free: 
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V.u> = 0. As we shall see this also caused problems, and again some methods had recourse 
to a Poisson equation (V 2 U = —V x u;). 

A more direct response to the question of the constraints is to “build them” into the dis- 
cretization. For instance the Contour Dynamics method (Zabusky, Hughes, and Roberts, 
1979) can be regarded as based on computational elements that “carry” Vu>, the gradient 
of the ‘2-D (scalar) vorticity, but the elements are automatically arranged on chains so that 
the vector field that purports to be is indeed curl-free. The equivalent for the vortic- 
ity is the filament method, in 3-D. Such “chained” methods are inconvenient in terms of 
bookkeeping (note how the Contour Dynamics method was never applied to more than a 
few chains), but are more attractive intuitively. This will be discussed again in the context 
of the 3-D methods. 
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2. POINT- VORTEX METHODS 


2.1. Basics 

Point-vortex studies have a long history. Interest in them is justified by the fact that 
the dynamics of point vortices provide exact albeit singular solutions of the incompressible 
Euler equations in 2-D. Assuming one can perform the time integration either analytically 
or numerically with high accuracy, one can generate nontrivial Euler solutions with a small 
computational cost and very little memory. 

Consider a collection of N vortices at positions Xj, j = 1 ,7V, with circulations Tj. The 
vorticity distribution is 

N 

u ( x ) = S(x-xj) (17) 

i = i 


where 8 is the 2-D Dirac distribution (6 is infinite at the origin, 0 elsewhere, and its integral 
is 1). Thus we have a collection of “spikes” of vorticity. Inserting (17) into (6) we get a 
sum: 


N 


U(x) = - e, x 


£ r *ifr 


X - Xk 


k= 1 


X k 


12 • 


(18) 


As an exercise one can verify that the velocity field given by (18) is divergence-free, is 
irrotational except at the points x*, and has a circulation T* around any contour enclosing 
only the k th vortex. 

The other ingredient, which expresses the dynamics, is Kelvin’s theorem. If we make 
each vortex follow the fluid and imagine a small contour around it we get the result that 
the circulation Tj of each vortex is conserved. Therefore, 



dxj 

~dt 


U(xj,f). 


(19) 

( 20 ) 


The formulation is now complete. The unknowns are Xj and Tj, j = 1 ,7V. The time 
evolution is given by (19) and (20), with (18) giving U in (20). Note that (18) is singular 
when Xj = x this contribution is simply omitted from the sum. In other words there is 
no self-induced velocity for the vortex. The circulations are constant, and the only task is 
to move the vortices with the fluid. 

An early example of point-vortex studies is Rosenhead’s (1931). As an exercise one can 
program the equations (using one’s favorite explicit time-integration scheme for now) for 
two vortices with various circulations of the same sign or of opposite signs. The vortices 
should orbit around their centroid (I^Xi + r 2 X 2 )/(Ei + T 2 ) unless their circulations are 
exactly opposite, in which case the pair moves in a straight line. 

Even though (17) looks like a simple linear superposition of solutions, it is not, because 
the vortices interact through (18). This is how the nonlinearity of the Euler equations 
enters the vortex equations. Note also that the equations (18,20) can be derived from the 
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Hamiltonian H = - I\,T* log j Xj -Xk\/4ir with TjXj and yj as the canonically conjugate 

coordinates and momenta (Aref, 1983). H corresponds to (16) with the self-energy of each 
vortex (which would be infinite) omitted; in that sense it is the “interaction energy”. 
Naturally, H is conserved by the dynamical equations. 

To be rigorous, the point-vortex method can only represent very special vorticity distri- 
butions, namely superpositions of Dirac distributions as in (17). These Dirac distributions 
are singular, and so is the velocity field they induce: it tends to oo like 1/jx — Xfc| as x 
approaches x*.. Such flows are not realistic in the sense that they have unbounded veloci- 
ties, infinite kinetic energy, etc. However, if one considers flows with small but nonsingular 
vortices the point-vortex method is still useful as we show below. 

An example of a realistic thin vortex is the “spreading line vortex” which is a classical 
exact solution of the 2-D Navier-Stokes equations (Batchelor, 1967, p. 204). The vorticity 
is given by 


u>(x, t) = exp(— - — ), 

v ’ 4 nvt v A ' A ’ 


4 vt 


( 21 ) 


the vortex being centered at the origin (x_, = 0). As the product vt tends to 0, (21) 
approaches the Dirac distribution (17). The velocity induced by this vorticity is 




( 22 ) 


As vt tends to 0 this velocity tends, point by point, to the one given by (18). Thus if 
the distance between vortices is large compared with their radius (y/vi) the equations 
(18, 19, 20) will not be affected. More generally, if we have small patches of vorticity 
(not necessarily of the type of (21)) of extent e much smaller than the distance d between 
patches the interaction of the vortices through (18) will not be affected. The patches also 
have internal dynamics; they deform, and spread if viscosity is present. However it can 
be shown that the centroid of the vorticity in the patch is not affected by the internal 
dynamics (convective or viscous) (see (14)). Therefore if x_, represents the centroid of the 
patch, the self-induced velocity is still exactly zero. 

The conclusion is that as the ratio t/d tends to zero the motion of the patches tends 
to that of point vortices. This was shown rigorously by Marchioro and Pulvirenti (1983) 
for short times, and it is reasonable to expect that it holds for all times, except maybe for 
some very special initial conditions. For the gravitational law and well- separated planets 
the analog is the familiar point-mass approximation. Saffman and Meiron (1986) present 
a different argument, showing that point vortices provide a weak solution of the vorticity 
equation (7) with v — 0 (although according to C. Greengard and Thomann (1988) this 
is controversial, and depends on one’s interpretation of how the Dirac distribution acts 
on a singular function such as (18)). Thus the point- vortex method is not strictly limited 
to vorticity fields of the type of (17); it is also a good approximation for the motion of 
small, well-separated patches of vorticity. This is still too restrictive for most practical 
applications. Consequently the recent work with point vortices tends to focus on math- 
ematical issues like the Hamiltonian character of the equations, the existence of steady 
configurations, integrable and chaotic motion (Aref, 1983), the statistical properties, and 
mixing of the fluid surrounding the vortices. A more elaborate model for well- separated 
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patches is the “moment model” of Melander, Zabusky, and Styczek (1986), which includes 
an approximation of the inviscid internal dynamics. 

2.2. Difficulties with Vortex Sheets* 

Given sufficient computing power it is natural to use more and more vortices to obtain a 
finer and finer description of a given vorticity distribution, just like one refines a grid. This 
idea was first applied to vortex sheets. Such flows are not as strongly singular as point 
vortices: the vorticity is still singular but only like a 1-D Dirac distribution, the velocity is 
discontinuous but finite. In fact Rosenhead’s work was aimed at a vortex sheet, but he did 
the computations by hand and could not use a large number of vortices. When this became 
possible, the outcome was a disappointment: smooth sheets could not be maintained, and 
the vortices rapidly started moving chaotically (Birkhoff, 1962). The time integration was 
accurate, and the situation did not improve as more vortices were used (Moore, 1979). 
Ostensibly, the method was not converging. 

The consensus now is that the wrinkling of the sheets discretized with point vortices was 
not just a numerical phenomenon, because the vortex-sheet problem is not well posed: in 
general the exact solution does not remain smooth even if the initial condition is. This 
is so because waves of all scales are unstable, and the shortest waves grow the fastest. 
Krasny (1986) and other researchers have shown that a periodic vortex sheet disturbed 
by a single sine wave remains smooth for a finite time, and then develops a singularity 
(infinite curvature). Krasny explains that if the initial condition is analytic the amplitude 
of the short waves is exponentially small, so that even though they grow very fast they 
do not contaminate the whole solution until some finite critical time t c . Krasny showed 
that the point- vortex method gave a good solution up to time t c , but only if round-off 
errors were reduced to a low enough level. He suggested using higher-precision arithmetic, 
and also proposed a filtering technique which is very powerful, but should be used with 
exquisite care. 


2.3. Application to Smooth Euler Solutions* 

The point-vortex method is not considered a viable numerical method to converge to 
smooth solutions of the 2-D Euler equations. However, the reason is not that it would fail 
in a blatant way, as it seemed to do with vortex sheets. The main reason the point-vortex 
method is out of favor is that convergence proofs have been given for vortex-blob methods 
(§3) and demand that the blobs overlap, and overlap more and more as they become 
denser. This is incompatible with point vortices. There is the question of what measure 
of the error one is using. If one takes the velocity or vorticity fields, there is no way the 
point-vortex solution could converge in any of the familiar norms, e.g., since it is not 
even bounded. On the other hand if one only considers the trajectory of the vortices and 
their initial arrangement is regular, they behave well at least for some time, which implies 
that the global measures like the moments of the vorticity and the Hamiltonian (13-16) 
behave well too. It may be that the trajectories (and even the velocity field if one takes 
an “exotic” norm like a Sobolev norm with a negative index) converge to the exact ones 
at least up to some critical time. This is just a conjecture. 
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3. VORTEX-BLOB METHODS 


3.1. Motivation and Basics 

Initially a strong motivation to switch from point vortices to nonsingular (or less- 
singular) elements than was the failure of the method on vortex sheets (Chorin and 
Bernard, 1973, Kuwahara and Takami, 1973). Now, the emphasis is on smooth Euler 
solutions, for which the blob methods have been shown to converge, even with high orders 
of accuracy (Hald, 1979, Beale and Majda, 1982). We present the method in 2-D. In a 
vortex-blob method the Dirac distributions 8 in (17) are replaced by bounded functions 

K: 

N 

w(x) = S T > ( 23 ) 

i=i 

The 8 a function has an integral of 1, and the length scale <r (the “core radius”) is a 
measure of the spread of 8 a . In almost all cases, 8 a is radially symmetric and its shape is 
independent of <r, i.e., 

M* - *<) = ^ (24) 

for some nondimensional function /, usually bell-shaped, and satisfying 

2rr f rf(r)dr = 1. (25) 

Jo 

One now has a “realistic” velocity field, bounded and with locally finite kinetic energy. 
The function / is small away from the origin, and as a — > 0, the blobs resemble point 
vortices. The Gaussian, as in (21), is a possible choice but may not be the most judicious 
in terms of computational speed. Equation (18) becomes 


U(x) 


N 


2tt Cz X 




x-Xfc 


k—i 


x* 


( 26 ) 


where F is defined by 

F(r) = 2tt f r'f(r')dr'. (27) 

Jo 

Figure 2 gives an example of the velocity given by (26, 27), compared with the point- vortex 
velocity (18). They agree except for \x\/(T of the order of 1 or less, and the blob-induced 
velocity is smooth at the origin. Finally the stream function is given by 

*w=sE r ‘^^i < 28 ) 

k = 1 

with F defined by 

Hr) = f dr'. (29) 

Jo r 
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Figure 2. Azimuthal velocity induced by vortices. — vortex blob; - - - point 
vortex 


Equation (19) is unchanged; circulation is conserved. In most methods (20) is also 
unchanged, i.e., the blob moves with the velocity at its center Xj. Leonard (1980) discusses 
the use of the velocity averaged over the blob (weighted by 6 a ) instead of the velocity at 
the center. In that case, (27) does not hold and the relation between / and F is slightly 
more complex; the alteration of the velocity is 0(cr 2 ). Averaging the velocity is sensible 
since the most meaningful interpretation of Xj is as the centroid of a small vortex patch 
(§2.1). Leonard also shows that the kinetic energy of the flow (16) is conserved only with 
the weighted scheme. On the other hand, this energy strongly depends on the function 
/, which has little physical meaning, and given a function F one can always define a 
Hamiltonian which approximates the kinetic energy, and will be conserved. The weighted 
scheme has not received much use. 

One could consider using a different core radius <Tj for each vortex; however this is ill- 
advised because it prevents the method from conserving even the most basic invariants, e.g., 
the linear momentum (14) (Leonard, 1980). One can of course symmetrize the interaction 
by taking an average of <Tj and <r*., which restores the conservation of (14), but this is 
arbitrary and seems inconsistent with the concept that the exact value of <Tj is important. A 
much better response to the objection that (14) is not conserved would be to use Leonard’s 
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weighted scheme. Still, the consensus is that the calculations are, and should be, quite 
insensitive to the function / and the value of <r. The basic method defined by (19, 20, 26) 
with uniform a (independent of j or t) is by far the most used. 


3.2. Convergence in Theory 

The accuracy and convergence have been discussed by Hald (1979), Leonard (1980), 
and Beale and Majda (1982), among others. Leonard focused on the consistency of the 
method whereas Hald, and Beale and Majda obtained complete convergence proofs includ- 
ing stability. Both approaches produced “moment conditions” which indicate the order of 
accuracy of the method according to whether some integrals of the function / equal 0. 

Beale and Majda split the error in the velocity field at some time t into two components, 
one due to the inexact position of the vortex centers (compared with the motion of the 
points in the exact solution), and the other due to the discretization of the velocity by 
(26). The two components naturally interact as t increases. Leonard examined the local 
error due to the fact that the whole blob is made to move at the same velocity, whereas in 
reality the fluid region that coincides with it is strained. He did not emphasize the error 
due to the initial discretization. 

The mathematical papers all show that convergence is obtained when N — + oo, h — * 0, 
and a — + 0, but a not as fast as h (i.e., cr/h — > oo). This last condition was unexpected 
(Leonard, 1980); it means that the vortices must become closely spaced, but also overlap 
more and more for the method to converge. Empirical tests by Nakamura et al. (1982) 
and Perlman (1985) fully confirm the need for the increased overlap. Hald and Beale and 
Majda suggest simple laws like <r = Ch q with C some constant and the power q less than 1. 
The studies by Hald, Leonard, and Beale and Majda all show that the error is proportional 
to some power of <r, not h. For instance with a Gaussian core function / the error is of 
order a 2 . 

The initialization of the vortices is more involved than one might expect , even in the 
simple case of an unbounded domain. Suppose one has an initial vorticity field w 0 . One 
lays a grid of cells of size h over the domain. One procedure is to place a vortex at the 
center of each cell, and to give it the circulation within that cell (the integral of u> 0 over the 
cell). The vorticity centroid would then be preferable to the geometric center of the cell. 
Another is to give it. the value of w 0 at the center, multiplied by h 2 . The first procedure 
is more intuitive and is preferable if u>o is not smooth; the second one is formally more 
accurate and is recommended if one is seeking better than second-order accuracy and u> 0 
is very smooth. A third procedure is to require that the vorticity given by (23) equal 
wo, either at the grid points or at a denser set of points with a least-squares algorithm 
(Nakamura et al., 1982). This results in a matrix equation for the Tj’s. This idea is 
attractive, but the implementation can be awkward, and the matrix is very ill-conditioned 
if (r/h is large. 


3.3. Convergence in Practice 

Convergence has been tested and observed only for unbounded flows (wall-bounded 
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flows are discussed in §0.0). Many papers contain some useful information about the 
behavior of the method (Nakamura et al ., 1982; Beale and Majda, 1985), but Perlman’s 
1985 study is the most complete. She chose some simple, smooth, radially symmetric 
Euler solutions and systematically computed the errors in actual calculations, including 
the various components introduced by Beale and Majda (1982). She found that for smooth 
u>o the method did converge, with the expected order of accuracy. For less-smooth u>o 
the liigher-order cores (ones which satisfy more moment conditions) did not provide any 
advantage. Values of the power q close to 1 (little overlap of the cores) made the calculation 
lose its accuracy in a shorter time. The disorganization of the vortices, compared with 
their regular initial arrangement, was largely responsible for the loss of accuracy. 

The question of disorganization is very relevant, because in practice the vortices always 
look disorganized, see figure 1. In all flows of interest, the product of the time and the 
velocity gradients (strain rate, vorticity) is large, so that the fluid gets highly distorted. 
The advantage of the vortex method over the few grid-based Lagrangian methods that have 
been tried is precisely that it can tolerate this distortion. However, its order of accuracy 
seems to drop to lower values. Quite probably, most vortex calculations produce results 
that are quantitatively accurate (especially the global quantities, e.g., forces and moments) 
but are far out of the range where the theoretical error estimates apply (§1-3, 9.6). This 
is at the same time a weakness, and a strength. 
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4. THREE-DIMENSIONAL FLOWS 


4-1. Essential Differences with Two-Dimensional Flows * 

There are two differences: the constraint on the vorticity to be divergence-free is not 
trivial, and infinitely thin vortices are not a viable approximation. The constraints were 
discussed in §l-4- In 3-D an arbitrary collection of Dirac functions as in (17), or of their 
smoothed counterparts as in (23), does not satisfy V.u> = 0 and therefore cannot represent 
a vorticity field. This difficulty has been addressed in two ways. One is to chain the vortex 
elements together to form filaments or nets, which automatically imposes the constraint. 
The other is to ignore the constraint at the level of the discretization and to rely on the 
accuracy of the method: as the method converges the constraint is better and better 
satisfied, as is the other equation (4). The first approach produced the filament and 
segment methods, while the second produced the “monopole” methods. 

The second difference is also a major difficulty. The velocity of a thin vortex filament of 
circulation T, radius of curvature R, core radius a, and binormal vector b is approximately 
given by 

u *sW!) b < 30 > 

(see Batchelor, 1967, p. 523). Thus filaments with curvature have a self-induction, and 
moreover it is not well behaved in the limit of an infinitely thin filament ( a/R —* 0). The 
next term in the approximation, after (30), strongly depends on the shape of the vorticity 
distribution, e.g., a “top hat” or a Gaussian as in viscous flows (21). This behavior of 
the self-induced velocity has motivated a number of studies based on the Local Induction 
Approximation (LIA) (Hama, 1962). Consider a smooth vortex with circulation T and 
global length scale L, and assume that a/L — ► 0. If one integrates from time 0 to T and 
the product T T\og(L/a)/L 2 is kept constant, the motion of the vortex tends to a finite 
limit as a/L — > 0, given by (30). Observe that (30) depends only on the local characteristics 
of the vortex. 

The LIA is reminiscent of the 2-D point-vortex approximation since it describes the 
motion in the limit of very thin vortices, but is in fact almost the opposite since the 
internal dynamics, which were negligible in 2-D, now control everything. Also, two distinct 
filaments do not interact at all, and different parts of the same filament interact weakly 
through the derivatives that enter the definition of the tangent, the binormal, and the 
radius of curvature. The LIA became especially popular when Hasimoto (1972) showed 
that it could support solitons. However in practice it is not easy to justify the value of a or 
to decide whether a should vary in space or time. The restriction a « L is also unrealistic 
in any situation with viscous or turbulent stresses. At this point this approximation has 
given excellent results in a few specific cases (Leonard, 1985), but it is not regarded as a 
numerical method adapted to the solution of practical problems. 


4.2. Filament Methods * 

Filament methods have been introduced and refined by Leonard and his co-workers over 
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the years (Leonard, 1985). Since vorticity is naturally arranged on closed lines Leonard 
chooses a discretization based on vortex filaments each with a constant circulation T j . The 
filaments then induce a velocity field on each other and on themselves, and their motion is 
computed. The self-induced velocity requires special care because the Biot-Savart kernel, 
equation (3), is strongly singular and many different regularizations have been tried. The 
method is in a state of flux, because the regularization that has been used up to now is 
accurate for long waves (it agrees with (30)) but is not satisfactory for short waves on 
the filament, often resulting in a spurious instability (Leonard, 1985). Recent work by 
Winckelmans and Leonard (1987) focuses on a solution to this problem which actually 
borrows a term from the LIA and adds it to the usual regularized contribution. 

Leonard’s method is by far the most used in 3-D (Leonard, 1980, 1985, Ashurst, 1983, 
Ashurst and Meiburg, 1985, Nakamura, Leonard, and Spalart, 1986). See also the closely- 
related 3-D vortex-in- cell method of Couet, Buneman, and Leonard (1981). With proper 
care this method can be applied to a variety of vorticity-dominated flows. Some involve 
isolated filaments, in which case the core parameters have a strong influence on the solution; 
others involve “bundles” of numerical filaments to describe a physical vortex, which makes 
the calculation much more expensive but less sensitive to the parameters. A systematic 
application to flows past solid bodies, with boundary layers, separation, and so on, still 
has to be made. 


4-3. Segment Methods * 

Chorin (1980, 1982) proposed a method in which chains of vortex elements can also be 
identified, but they are treated totally separately in terms of the velocity they induce, and 
the simplest and least computationally expensive regularization of the kernel is used (the 
velocity is corrected within a sphere centered on the mid-point of the segment). Chorin 
does not require a single chain to behave like a physical vortex, and relies on bundles or 
clouds of elements to obtain the correct behavior (e.g., (30)). Compared with filament 
methods this method is somewhat cheaper computationally and offers easier bookkeeping. 
In fact, one could generalize the pattern of elements from a bundle of lines to a “net”. One 
would keep track of a set of nodes, which get moved in time following the fluid. There 
would also be a set of segments, each with a circulation, a starting node and an end node. 
Any number of segments could start from, or end on, a single node provided that their 
circulations (taken with the appropriate sign) add up to zero. With such a method it 
would be much easier to merge nodes to control their number if needed (§9.3). 

Chorin (1982) used his method to compute the evolution of a disturbed vortex in a 
periodic box, with emphasis on the question of singularities in the 3-D Euler equations 
and the fractal character of the support of vorticity in turbulent flow. He also studied 
boundary-layer transition (Chorin, 1980). 


4-4- Monopole Methods * 

A number of 3-D methods have been proposed in which the computational elements 
are totally disconnected instead of being arranged in filaments or nets (Rehbach, 1977; 
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Saffman, 1980; Beale and Majda, 1982; Novikov, 1983; C. Greengard, 1984; Cantaloube 
and Huberson, 1984; Anderson and C. Greengard, 1985; Mosher, 1985). The elements 
are sometimes called “vortex monopoles” or “vortons.” In such methods the rotation and 
stretching of vorticity by the u;.VU term is explicitly computed instead of being a by- 
product of the relative motion of the points that describe the filament. The methods differ 
primarily in the way this term is computed. 

The attractive features of such methods are: simplicity and computational convenience, 
possibly easier mathematical analysis, and the freedom for the elements to rearrange them- 
selves so that the topology of the vortex lines changes (the “reconnection process”). The 
objection was raised early on (by the author and certainly by many others) that in such a 
method the vorticity field carried by the computational elements may not be divergence- 
free, which is unnatural. 

Suppose one has a set of elements which carry vorticity in 3-D space. Call u> a the vector 
field formed by adding the vorticity of all the elements, as in (17) or (23) but with vector 
values. The subscript “a” stands for “apparent”, as will become clear. The procedure is 
to apply the Biot-Savart law (3) to to obtain the velocity field U. Let us define the 
“true” vorticity field by w ( = V x U. Clearly the divergence of u> t is 0, since it is the 
curl of some vector field. If the divergence of the apparent vorticity is nonzero, it cannot 
be equal to the true vorticity: o> a ^ u;*. More precisely, u> t is the projection of u; a onto the 
space of divergence- free vector fields (this is so because the Biot-Savart law actually solves 
the Poisson equation V 2 U = —V x u>„, so that the only link between u> a and is that 
they have the same curl). 

Graphically, what this means is that when the “vortex lines” of u a terminate within 
the domain (i.e., V.a>„ ^ 0) the vortex lines of u? t correct this by connecting to another 
element, or the opposite end of the same element, in as smooth a fashion as possible. 
Thus the small region in which u> a / 0 is surrounded by a much larger “halo” in which 
u> t 7^ 0. This is disturbing, because one applies the vorticity transport equation (including 
the vortex-stretching and rotation term, etc.) to u> a , not to u) t . Thus if there is rotation 
locally where uj a is located, the whole halo of w t rotates too, which is just not correct. 

If the visualization of a monopole calculation shows that the elements are scrambled 
and some of the vortex lines of u> a terminate within the domain, then w a and must be 
strongly different, the calculation is unreliable, and it should be terminated. Monopole 
methods are unlikely to be “robust” in the sense of allowing safe integration to large times. 
Rehbach (1977) reports serious problems far downstream of his wing and had to artificially 
reduce the strength of the vortex elements in that region to obtain a stable solution. 

Saffman and Meiron’s (1986) study is quite critical of vorton methods. They showed that, 
in general, a system of vortons is not a weak solution of (4) and does not even conserve 
the integrated vorticity or momentum, given by (9) and (10) (see also Leonard, 1985); 
they also comment on the central role of the nonzero divergence of u?„. Nonconservation of 
integrated vorticity is not too serious, because it is the integral of u? a which is not conserved 
( u> t always integrates to 0, assuming fast enough decay at large distances). On the other 
hand, nonconservation of momentum is a serious flaw (§1.5). 

The most popular application of monopole methods has been to the reconnection process 
(C. Greengard, 1984; Mosher, 1985), which filament methods are unable to model without 
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outside intervention. Rehbach (1977) also studied the flow past sharp-edged 3-D wings, 
and Cantaloube and Huberson (1984) studied propellers and rotors. Only Greengard 
had a mathematical convergence proof for his method and attempted (with some degree 
of success) to demonstrate the convergence in actual calculations. The 3-D problem is 
probably the most open and challenging in the field of vortex methods. 


5. EXTENSION TO VISCOUS AND COMPRESSIBLE FLOWS 


5.1. Viscous Flows 

The original Navier-Stokes equations express local interactions between neighboring 
points, through the spatial derivatives. In contrast, vortices interact at any distance 
through the Biot-Savart law and do not recognize their neighbors. The unfortunate con- 
sequence is that while we were able to formulate one set of equations, the Euler equations, 
very elegantly the necessary spatial derivatives are not available to add any other terms, 
for instance the viscous terms, in a straightforward manner. One needs to alter (19), (20), 
or (23). 

One obvious idea is to use the exact viscous solution (21) to make vortex blobs, at least 
in 2-D. The blobs then spread in time like \fvt. Equations (19) and (20) are retained. This 
method is attractive, but inconsistent. In a vortex- blob method, to obtain convergence 
of the transport process one needs to let the core size <7 tend to 0 (see §3.2). If <r is 
made proportional to \A 4, the user can’t control it any more, and the core-deformation 
error cannot be made to vanish. Thus the method “mimics” viscosity, but damages the 
inviscid process. C. Greengard (1985) actually proved that the method solves “the wrong 
equation”: the vorticity is “correctly diffused, but incorrectly convected, even in the limit 
of infinitely many vortices”. 

The most used treatment of viscosity is the random- walk algorithm (Chorin, 1973). 
Chorin adds a random component of appropriate magnitude (proportional to \/vAt) to 
the position of each vortex at each time step. In the average, with a large number of 
vortices, this introduces the desired diffusion. When solid bodies are present there is also a 
creation of new vortices along the wall at each step to account for the flux of vorticity out 
of the wall. This approach is very simple, consistent, and partial proofs of convergence are 
starting to appear (Hald, 1984, Goodman, 1987). On the other hand it seems inefficient 
to carefully track a vortex if it contains only such a small amount of information (since 
it takes many vortices for the laws of statistics to turn the random walk process into a 
meaningful macroscopic diffusion). Milinazzo and Saffman (1977) created a controversy 
by testing the method and showing very slow convergence to a simple exact solution; the 
error was of the order of iV -1 / 2 , with JV the number of vortices. There are some questions 
about which measure of the error is most meaningful with the random vortex method, 
but the consensus seems to be that slow convergence is observed notably near boundaries 
(Chorin, 1980), and the 0(JV -1//2 ) estimate was confirmed by Roberts (1985). 

Probably the most convincing application with random walk was the study by Ghonieni 
and Gagnon (1987). They succeeded in simulating the flow over a backward-facing step 
in the laminar range, up to Reynolds number 229, using about 1000 vortices. In this 
flow viscous stresses control the flow even away from the wall, so the agreement with 
experiments and other calculations is very encouraging. On the other hand, in the author’s 
estimation results of the same quality could be obtained with a grid-based method with 
much less computational expense, and the competitiveness of the random- walk method in 
terms of computational speed is not excellent. This led the present author to choose a less 
ambitious approach, wherein the viscous effects are totally neglected away from the wall, 
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and are included only in a relatively crude manner in the boundary layers. This approach is 
clearly not universal, but with careful implementation and high enough Reynolds numbers 
it has been justified by the results obtained on a number of cases. 

Two other possibilities for a viscous method should be mentioned. One could be called 
the “exploding vortex” method: at each step each vortex would be split into a few vortices 
with a spread again proportional to y/ uAt. This method leads to a rapid proliferation 
of vortices, and therefore is extremely expensive and impractical. Finally, there is the 
idea of transferring circulation between blobs at each time step in an appropriate manner 
to represent the diffusion (G. Winckelmans, personal communication, 1987, after work 
by Cottet and Raviart’s group in France). Thus, (19) is modified. There is no core 
spreading, so the convergence of the inviscid part of the dynamics is not jeopardized. 
Circulation is conserved, but the vorticity centroid apparently is not (one could modify (20) 
to compensate). The method needs to have some “dormant” vortices with zero circulation 
in any region that may become vortical at later times, and the vortices must be closely 
packed, only a distance of order y/uAt apart. In grid-based methods the spacing between 
points is typically of the order of yfvT (where T is a global time scale), which is much 
larger than y/uA t. 

The profusion of candidates for the viscous treatment (we even omitted a few that 
showed very little promise) reflects the intense need for such an addition, as well as the 
lack of a truly satisfactory solution to this date. 

5.2. Compressible Flows* 

Most flows of interest are at least slightly compressible, and the capability to treat such 
flows with a Vortex method would be very valuable. Furthermore one of the key ingredient s 
of the method, Kelvin’s theorem of conservation of circulation, still applies to barotropic 
compressible flows. Shocks and other sources of vorticity invalidate the theorem; this 
creation of vorticity could in principle be computed, but this has not been attempted yet. 
Let us restrict the discussion to shock-free flows. Then (19) and (20) apply. 

The equation that does not apply any more is (18). The velocity now also depends on 
the dilatation rate (u x + t'j,), which is nonzero. One possible approach is to include this 
quantity in the description, and have a number of sources in the flow just like there are 
vortices. The obstacle here is that the evolution equation for the dilatation rate is far from 
being as simple as that for vorticity. Furthermore the dilatation rate, unlike the vorticity, 
is not normally confined to a small region. Thus one of the main advantages of the vortex 
method is lost. These two obstacles seem to have prevented any sustained attempt using 
this approach. The author is only aware of some work along these lines by Professor J. C. 
Wu, using a finite- difference method (personal communication, 1983). 

Another approach is to consider flows at low Mach numbers and to approximate the 
dilatation rate. One would then obtain an expansion in terms of the Mach number, similar 
to the well-known and successful Prandtl-Glauert approximation. This was attempted 
by Deffenbaugh and Shivananda (1980). They obtained a Poisson equation for the first 
correction, of the order the Mach number squared, and used a finite-difference method 
to solve it. It is not clear to the author that the procedure was consistent, because time 
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derivatives were found both on the left- and right-hand side of the equations. Indeed the 
numerical solutions were unstable in time, and Deffenbaugh and Shivananda had to delete 
some of the time-dependent terms to obtain bounded solutions. They conjectured that a 
more accurate scheme would alleviate the problem. To our knowledge, this approach has 
not been followed upon. 

For future work, it seems that the objections to the first approach are deeper than those 
to the second one. One should be able to formulate a correction at low Mach number. 
However it is not clear that this can be done without the use of spatial derivatives as in a 
Poisson equation; these can be handled on a grid, but losing the grid-free character of the 
method is a heavy price to pay. The best by far would be to find a grid-free strategy for 
the compressibility correction. 
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6. INVISCID BOUNDARIES 


6.1. Use of Boundary Elements 

The traditional way of implementing the inviscid boundary condition (i.e., zero normal 
velocity at a surface) is by using image vortices. However this is possible only when the 
shape of the surface is very simple (e.g., straight line or circle) or when one has an analytical 
conformal mapping from the desired shape to one of the simple shapes. This has led to 
numerous studies of the flow past circles, ellipses, and Joukovsky airfoils, for instance. 
Another serious drawback of images in practice is that a vortex close to the wall is also 
close to its own image. A strong interaction can develop, with the vortex and the image 
forming a pair and traveling along the wall at a high velocity and in the direction opposite 
to that of the boundary layer. 

In practice the method of boundary elements, in which vortices, sources, or doublets of 
appropriate strength are placed along the boundary to ensure that the boundary condition 
is satisfied, is much more flexible. In other contexts it is called the panel method. At the 
expense of solving a linear system, one can treat arbitrary shapes. The matrix of the linear 
system is unfortunately full, but in a time- dependent computation it needs to be inverted 
only once so that the cost at each step is of order M 2 and not M 3 (M being the number 
of boundary elements). 

When the boundary-element method is used as part of a vortex algorithm, as opposed 
to a potential-flow solver, it is convenient to use vortices for the boundary elements (Lewis, 
1981). Furthermore in some cases using vortices is much more natural (for instance the 
splitter plate for a mixing layer carries a velocity jump and therefore circulation). Let us 
restrict our attention to 2-D for now. Let. the points 1 < j < M , describe the boundary. 
The bound vortices are placed at these points. One could construct an algorithm to cancel 
the normal velocity at or near the points Xj (for instance at (x.j + Xj+])/2). However, there 
are large local velocities and depending on the exact choice of the test point and other 
parameters there is a danger of having fluid “leak” through the boundary in the average. It 
is much preferable to use the stream function: the condition is then V’(Xj+i ) = V’(Xj)- Thus 
we are imposing zero mass flux between Xj and Xj+j. This method is more robust, and 
also more convenient when one has several surfaces (e.g., the walls of a wind-tunnel) and 
wishes to impose the mass flux between them: the condition becomes ip(xj + j ) = V’(x? ) - Q 
if Xj is on one wall, Xj+y is on the other and Q is the mass flux. The stream-function 
method results in a matrix which is well structured; the largest elements are on the main 
diagonal and the one just below it, and Gauss elimination without pivoting is sufficient to 
solve the system. 

Typically the stream function results from the linear combination of three components. 
If there is flow at infinity this corresponds to a stream function ^> 00 ( x ) = (x x Upoj.ej 
where U,*, is the uniform freestream velocity. Note that a slight generalization is possible 
in 2-D: if the freestream flow has uniform vorticity the vortex method can still be applied 
to the deviation from the freestream vorticity (i.e., equation (7) still holds, but equation (4) 
does not in 3-D). Then there is the stream function associated with the free vortices, V’/> 
which is given by (28). The initialization or generation of these vortices strongly depends 
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on the problem at hand. Finally there is the stream function of the bound vortices, 
which is given by (28) too. Thus the equations are 


V>b(Xj + |) - V’fc(Xj) = V’oo(Xj) - ^oo(Xj + l) + V\f(*j) - - Qr ( 31 ) 

The unknowns are the circulations T'j of the M bound vortices. After they have been 
computed, the velocity of the free vortices is computed and they are advanced for one time 
step according to (20). The displacement of the free vortices (and possibly the creation of 
new ones) alters V’/, so that (31) has to be solved again for the next step. 

The question arises immediately that (31) can be written only for j < M; there are 
only M — 1 equations, for the M unknowns 1^. This is the well-known paradox: the total 
circulation is undetermined. The solution is to add an M th equation that prescribes the 
sum of the circulations. Thus all the elements of the last row of the matrix are 1. Typically, 
one forces the total circulation in the plane to be zero, so that the M th equation is: 

M N 

E r ; = -E r < < 32 > 

j = 1 i= 1 

(recall that N is the number of free vortices). However this condition is not unique; 
sometimes symmetry or other considerations are applied instead. See the discussion of 
(38) in §7.2. 

For the extension to 3-D, one will not have a stream function. One can of course revert 
to a condition on the normal velocity itself, but imposing the condition in an integral sense 
as in (31 ) is preferable. Suppose one has cells on the surface (e.g., triangles); the border of 
each cell defines a contour. One could use a vector potential and require that its circulation 
around each contour is zero. This is equivalent to requiring that the flux of the velocity 
field through each cell is zero. Since mass conservation is built into the discretization, it 
seems consistent to use for the boundary an algorithm that exactly conserves mass. 

An application of the vortex method with inviscid boundaries is given in §6.3. Other 
applications are nozzles or wind tunnels. Spalart (1984) assessed tunnel blockage effects 
by imposing the inviscid condition on the tunnel walls (on which separation did not occur) 
and viscous conditions on the body. One can also represent a lifting surface as an inviscid 
boundary (Katz, 1981), but then the circulation in the lower and upper boundary layers 
cannot be distinguished, so that one cannot predict viscous effects such as separation; one 
also needs to impose the Kutta condition. 


6.2. Remark on the Choice of Two Numerical Parameters * 

In general, there are no explicit equations linking the different length parameters, so 
that the user is left with his/her intuition to choose values. However in the simple case 
of an inviscid boundary one can use the following simple argument to relate the spacing 
between the points Xj and Xj +1 and the core radius a. We idealize the boundary as a 
straight line on the x axis with equally distributed points Xj at intervals A. Let U denote 
the velocity jump across the sheet. The ideal stream function is V’i = — U\y\/2; the flow is 
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uniform at ±f//2 above and below the sheet. With the discretization by vortices, the flow 
is nearly uniform away from the sheet by forms the usual “cat’s eyes” around the vortices. 
The stream function discretized using vortex blobs with the core function F = r 2 /(r 2 -f 1) 
(see (46)), is given by 
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The circulation of each vortex is — U A, and the denominator inside the logarithm represents 
an additive constant that was set to ensure that at x = Xj = j Ae x , = ipi = 0 (recall that 
the stream function will be sampled at the points Xj to write (31)). One can verify that 
for large |y| the stream function tpd given by (33) becomes arbitrarily close to — U\y\/2 + C , 
where C is a constant. The mass flux between the test point Xj and a point far up will 
be correct only if C = 0. This mass flux is an important quantity when the walls are the 
walls of a tunnel or nozzle, especially of variable section. It is less important in external 
flows. A numerical summation of the series (33) showed that C = 0 is obtained when <r/A 
is equal to 0.153. The optimum value of <t / A naturally depends on the core profile used 
(the F function). In practice, the value 0.153 can rarely be exactly respected since the 
spacing of the points is generally uneven, but it is a useful guide. 


6.3. Example: a Curved Mixing Layer 

Mixing layers, either temporally- or spatially-developing, have been a frequent subject 
for vortex methods (Ashurst, 1979; Aref and Siggia, 1980; Ashurst and Meiburg, 1985). 
A short, unpublished study of curved, spatially-developing mixing layers will be described 
here primarily as an example of how realistic boundary conditions can be applied. The 
study was suggested by Dr. N. Mansour (NASA Ames Research Center) after some exper- 
imental studies showed a strong effect of curvature on the structure of mixing layers (with 
the faster stream inside, the layer was much more disorganized and three-dimensional). 

Although this was not the only way to proceed (see the treatment of the flat mixing 
layer using two semi-infinite vortex sheets in Leonard, 1980) it was decided to use curved 
walls as guides for the flow. Let R be the radius of the mixing layer, Rj the radius of 
the inner wall and R 2 that of the outer wall. The values chosen satisfied R\ / R = 0.9 and 
R2/R = 1.1, and the “test section” covered one radian (see figure 3). Boundary conditions 
were enforced on the two walls, on an inflow boundary and on a section of the splitter plate 
using bound vortices. On the walls and plate, the inviscid no-penetration condition was 
applied (most mixing-layer simulations do not enforce any condition on the plate; they just 
place a row of vortices of equal circulation on it). Along the inflow boundary, the velocity 
normal to the boundary (uq) was imposed and equal to the velocity of the irrotational flow 
(i.e., proportional to 1/r). This is an example of (31) with Qj ^ 0. It can be interpreted 
as a porous boundary with imposed blowing. The system was closed by requiring that 
the circulation of the bound and the free vortices add up to zero, see (32). This value 
was arbitrary but ensured that the flow outside the walls was relatively quiescent, so that 
the velocity field was smooth at the inflow and outflow boundaries. Provided that the 
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Figure 3. Computation of curved mixing layers, a) Fast fluid outside, Vi — W\\ 
b) fast fluid inside, Ui = V\ /2. o free vortices; + vortices bound to 
walls; A vortices bound to splitter plate; □ vortices bound to inflow 
boundary 





test section is long enough, the manner in which the linear system is closed (the M th 
equation) is unimportant, as are the inflow and outflow conditions; here we are invoking 
Saint- Venant’s principle. 

No condition was needed at the outflow. The trailing edge of the plate gives an example 
of a smooth transition from the bound vortices to the free vortices; the two families of 
vortices have the same spacing, core size, and essentially the same circulation (this assumes 
the proper relationship between the mean velocity, the time step, and the spacing). At 
every time step the last bound vortex becomes a free vortex, keeping its circulation. The 
vortex that is farthest downstream is deleted. Thus the calculation can be run as long as 
desired. 

In parallel with the numerical study, the linear stability of curved vortex sheets was 
investigated. By generalizing the well-known derivation for flat sheets (Batchelor, 1967, p. 
511) it was found that the growth rate a for a wave with wavenumber a is given by: 

a ! = ~(V t -U,) 2 + ] ^Ui-Ul) (34) 

where U\ is the velocity of the inside stream and t/ 2 is the velocity of the outside stream. 
Note that this formula reduces to Batchelor’s result (a = a\U\ — 2 1/2) as aR — > 00 , that 

is, as the radius of curvature becomes much larger than the wavelength. 

Equation (34) shows that the mixing layer is more unstable if the inner stream is faster 
(tf? > f/|). However the relative difference in growth rates is smallest for the most 
unstable waves (largest |a|) so this effect is not significant. Figure 3 shows mixing layers 
with velocity ratios U 2 /U 1 of 2 and 1/2 respectively, and confirms that even a significant 
curvature (compared with the size of the largest eddies) has little effect on 2-D mixing 
layers. This strongly suggests that the third dimension is necessary to explain the effect 
observed in experiments. 

Another point of interest is that the boundary condition on the splitter plate did not 
seem to make a noticeable difference in the flow pattern. There was a possibility that the 
pressure disturbances from downstream would influence the circulation of the new vortex 
that is created at each step, thus creating a feedback mechanism with the sharp trailing 
edge as an amplifier. However, in the calculations this effect was very weak, and the flow 
did not differ noticeably when the boundary condition was enforced using (31) and when 
the vortices bound to the plate were given uniform circulations. 
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7. VISCOUS BOUNDARIES 


7.1. Application of the Biot-Savart Law Inside a Solid Body 

In §6 we showed how the inviscid, no-penetration condition can be conveniently enforced 
using boundary elements, preferably vortices. We are now interested in satisfying also the 
no-slip condition. A priori this is more difficult, but we have a very fortunate result that 
shows that, for external flows, this can still be achieved with a single set of boundary 
elements (Wu and Sankar, 1980). Thus, calculations with both the no-penetration and 
no-slip conditions are almost as simple as those with only the former. The key idea is to 
extend the velocity field and the stream function into the region occupied by the body. 

Consider a solid body occupying the bounded region D, with boundary dD. It may 
be multiply-connected, made up of sub-regions Di. For now we restrict ourselves to a 
fixed body and 2-D. Suppose we obtained a vorticity field such that the stream function it 
induces through the Biot-Savart law (8) satisfies the no- penetration condition, drl’/ds = 0 
on dD. Recall that (8) represents the Green’s function for the whole plane. We may have 
achieved this using boundary elements; in any case there are only vorticity elements (so 
that ip is defined and single-valued), and none of them are inside D. The formula (8) can 
be applied for x inside D, and there V 2 ip = u> — 0. Therefore ip satisfies: 


V 2 V> = 0 in D, 

^ = 0 on dD. 

ds 


(35a) 

(356) 


This is Laplace’s equation with a Dirichlet boundary condition. It is a well-posed problem 
with a unique solution up to an additive constant in each subregion: ip = Ci in £>,. 
Therefore U, as given by (6), is identically zero in D. Also consider dip/dn , the normal 
derivative of ip as one approaches dD from the inside. Since ip is constant in Di , 


dip . 

— — = 0 on dD, 

dn 


(36) 


i.e., the tangential velocity just below the boundary is zero. This is a very useful result. 

Several remarks. First, one could have first imposed dip/dn — 0 and obtained dip /ds = 0 
as the by-product (Laplace’s equation with Neumann condition also has a unique solution 
up to an additive constant). In other words, if (35a) holds (35b) and (36) are equivalent. In 
the numerical method we elected to explicitly apply the no- penetration condition dip /ds — 
0 because it requires less smoothness of ip (observe that (31) is written directly in terms 
of ip, not in terms of derivatives) and also because it seemed more important to conserve 
mass exactly; the viscous effects, including the no-slip condition, are treated a little more 
loosely. Second, the solution to (35) is unique only if D is bounded and dD is closed, 
that is, for an external flow. Third, the argument on the normal derivative applies only if 
one can approach dD from the inside, that is, if D has a finite thickness. The argument 
totally breaks down for an infinitely thin body, e.g., a plate, and in the numerical context 
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Figure 4a. Flow past, a bluff body, just after impulsive start. — boundary; 
• vortices; — » velocity vectors 


one needs to be careful if the thickness of the body is not large enough compared with the 
spacing of the vortices along the wall. 

The behavior of the velocity field is illustrated in figure 4. A simple body shape was 
considered, just after an impulsive start (i.e., there are no free vortices yet). Thus we have 
irrotational flow around the body. Boundary vortices were used to enforce the drf/ds = 0 
condition, discretized as in (31). The velocity vectors computed using (26) are plotted 
at the nodes of a regular grid, including inside the body. In figure 4a one recognizes the 
usual potential-flow solution outside the body. Inside the body, the velocity vectors are 
very small, showing that the formal argument based on Laplace’s equation (35) does work 
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well numerically. This was not entirely obvious, because usually the function S a in (23) 
does not have bounded support, so that there is a little residual vorticity inside the body, 
and also because we forced i}> to take the same value (Cj) only at the discrete points Xj 
(see (31)). The natural smoothness of solutions to Laplace’s equation helps. In figure 4 
the boundary elements represent a vortex sheet, and the strength dT/ds of the sheet gives 
the velocity jump across the sheet. Therefore dr / ds gives the tangential velocity above the 
elements. Thus once the boundary vortices are known it is a very simple matter to obtain 
the limiting velocity of the potential flow at the wall, and the surface pressure through 
Bernoulli’s equation. Thus even for potential flow the method of boundary vortices has a 
distinct advantage over those based on sources or doublets. 

7.2. Creation of Vorticity 

The interpretation now differs for inviscid and viscous flow. In viscous flow the vortices 
along the boundary are not considered as mathematical tools to enforce a boundary con- 
dition, but as actual vorticity that is imparted to the fluid by the shear stress at the wall. 
Solid walls have long been interpreted as sources of vorticity (Lighthill, 1963). New vortices 
are created at every step and released into the fluid, just like at the trailing edge for the 
mixing layer (§0), but now all along the wall (Chorin, 1973). The circulation of the new 
vortices is computed as is they were bound vortices, but they are immediately considered 
as free vortices and allowed to move with the fluid. Symbolically, their contribution to 
is transferred from V’f> to V’/» on the right-hand side of (31), so that V’oo + V’/ alone satisfies 
the boundary condition. The vortices then move according to (20), by an amount that is 
O(At). This alters i/’oo + ipf by an amount O(At), so when (31) is solved again xfb and 
therefore the circulation of the new vortices are O(At). Thus, a small amount of vorticity 
is released at each step. The vortices are created not right on the boundary, but a small 
distance d 0 outside it. There are “wall points” Xj where (31) (with Qj = 0) is applied, 
and “creation points”, say x'j where the new vortices are placed. 

The wall points and creation points can be distinguished in figure 4b, which is a detail 
of figure 4a. Notice the smooth potential flow away from the wall, the velocity vectors 
becoming very small inside the body, and the transition region in which the direction of 
the velocity fluctuates. The microscopic description of the wall (on the scale of |xj + i -Xj|) 
is a little crude, but on a macroscopic level the method has the following advantage. The 
row of vortices in figure 4b (and in general with any attached boundary layer) forms a 
vortex sheet. The tangential velocity above the shear layer is U e , the velocity of the 
potential flow, and the velocity below it is zero. Therefore the vortices translate along the 
wall at a velocity U e / 2. This is known to be correct, in the sense that in a boundary layer 
the average velocity of the vorticity, defined by 

f°° 

U av —l u U) dy 

Jo 

is precisely equal to Z7 e /2 (just use u % — u y and an integration by parts). The flux of 
vorticity along the wall is an important quantity and is correct (as mentioned earlier, it 
can be very inaccurate when image vortices are used). This is another example of how 
small-scale errors can cancel each other so that integrated quantities (here, U av ) are correct. 
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Figure 4b. Flow past a bluff body, just after impulsive start (detail of Fig. 4a). 

-o- boundary and wall points x^; • vortices at creation points x^; 
— * velocity vectors 


As in the inviscid flows, we need an extra equation to close the system (31) for each 
body. This condition is that the circulations of the vortices created along each body at 
each time step add up to zero (Wu and Sankar, 1980): 

M 

^r; = « (38) 

J=1 

Note that this condition is rigorous, unlike (32) which looks the same but was arbitrary 
(actually, (38) was known and (32) was chosen by analogy). If one considers the whole flow, 
(38) implies that the total circulation is conserved, a well-known result, see (13). However, 
(38) is a stronger condition since it applies to each body separately. In summary, for each 
body described by M wall points there are M — 1 type-(31) equations and one of type (38). 
Note that the type-(31) equations are coupled between bodies, so that the whole matrix 
is full (except for the few type-(38) lines). 

In inviscid flow, u> can have a singularity of the order of a 1-D Dirac distribution, and ip 
is continuous but with a jump of its normal derivative across dD. In viscous flow, u> has 
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a jump across dD but is bounded, and both derivatives of ip are continuous. Therefore 
(36) applies on either side of dDi, i.e., the no-slip condition is satisfied. u> is bounded 
because with viscosity it cannot remain in a sheet; instead it diffuses into the fluid. The 
mechanism which prevents it from accumulating at the wall is viscous; thus the vortex- 
creation algorithm is consistent only with a viscous flow solver. In actual calculations 
numerical effects play a strong role in this process, see §7.^ and 8.2. 

Several slight generalizations are possible. If the body is translating, this can be taken 
into account using the Qj's in (31). If the body is rotating at an angular velocity f l the 
situation is a little more involved since there is vorticity = 2fl inside D. Again, one 
adjusts the Qj' s. The stream function induced by the inside vorticity can be computed 
with the Contour Dynamics formula (Zabusky, Hughes, and Roberts, 1979) and subtracted 
from the right-hand side of (31). This was done by Spalart and Leonard (1981) for the 
dynamic-stall study, see figure 7. On the other hand, if the body is deforming or several 
bodies are in relative motion there is an extra difficulty because the matrix involved in 
(31) becomes time- dependent, so that it has to be computed and inverted at each step. 
Finally, in 3-D the same argument as in (35) can be applied to the vector potential. Thus 
the extension to 3-D is possible, but it has not been made yet; we first need to obtain a 
robust method for unbounded 3-D flows. 


7.3. Choice of Numerical Parameters * 

In the creation algorithm just described, three length scales are involved near the wall. 
As before (§6.2), we have A, the spacing between the wall points, and <r, the core radius. 
Now there is also do, the distance from the wall to the creation points. The requirements 
conflict. d 0 should be large enough compared with a for the residual vorticity inside 
the body to be small, but it should also be small enough compared with A to maintain 
a strong coupling between the creation point and the wall point (so that the matrix is 
well-conditioned). Moreover, <r should be large enough compared with A for the cores to 
overlap in the s direction, so that the velocity field is as smooth as possible (figure 4b). 
The difficulty arises because the vortices are circular in a region, the boundary layer, where 
the length scales in the s and n direction are naturally disproportionate. This is partly 
what motivates separate treatments of the boundary layer (see §8). 

For viscous flows, a quantitative condition as in § 6.2 has not been obtained; ip 4 cannot be 
forced to give the right mass flux from Xj | to oo without making the vortex cores impinge 
on the boundary. In that sense, the actual boundary is blurred, somewhere between the 
wall points and the creation points. Figure 4b shows a satisfactory set of parameters: 
do = A a /2, & = A a /4, where A a is an average over the boundary: A = S/M where 5 is 
the length of the boundary. The need for a balance between A, d 0 and cr implies that the 
local value of A should not have very wide variations. It is helpful to cluster the points 
(thus reducing A) in the regions with high curvature and near sharp corners, but this 
clustering should be moderate. Typically, A should not vary by more than a factor of 2. 
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7.f. Kutta Condition 

The Kutta condition, like conformal mappings and image vortices, is a major tool of 
classical hydrodynamic theory, and it has been explicitly incorporated in many vortex 
methods (e.g., Stansby and Dixon, 1982). Actually, a generalization of the classical con- 
dition (which just provides the circulation for steady flow on an airfoil with a single sharp 
trailing edge) was needed since the flows have vorticity and are time dependent, and there 
may be several sharp edges, like on a square. However the idea is the same: the velocity 
at any sharp corner must be finite. Note that “Kutta conditions” have been applied to 
smooth surfaces at separation, which is probably an improper use of the word (in fact, one 
is just applying the no-slip condition). Note also that a complete viscous method does not 
need the Kutta condition, because this condition is simply the result of separation of the 
boundary layers, which an accurate viscous calculation can account for. 

It turns out that the method outlined in §7.2, with creation of vorticity and the no-slip 
condition satisfied all along dD , does a good job of reproducing this viscous behavior, even 
when the viscous terms are not included. For instance in figure 1 the flow near the various 
trailing edges and convex sharp corners is physically correct, although no special treatment 
or extra condition was applied there. When confronted with a convex sharp corner, the 
method naturally expels boundary-layer vorticity into the outer flow until a smooth flow 
without large velocities around the corner is established. The high curvature causes a rapid 
acceleration and deceleration of the boundary layer, which causes separation. Vorticity is 
released until the streamline originating at the corner is aligned with the upstream side of 
the corner. 

There are two major reasons for this “pseudo-viscous” behavior. The first is that, the 
vortices are created at a small but nonzero distance do from the wall so that whenever 
there is a deceleration of the boundary layer they acquire a small velocity component 
directed away from the wall, which is what initiates separation. The second is that time- 
integration errors, especially with large velocities and strong streamline curvature, move 
the vortices away from the wall. The analogies between integration errors and viscous 
effects are discussed in more detail in § 8.2 and 9.1. In any case this somewhat fortuitous 
behavior is consistent, and the flow near sharp edges has always been physically correct. 
Note that the shedding of vorticity can be only transient as on a starting airfoil without 
stall, or permanent as in figure 1 or in the case of flow past, a square. The separation from 
smooth surfaces is a much more delicate effect, and is discussed in detail in §8. 


7.5. Example: Starting Vortex on an Airfoil 

The flow development around an airfoil at incidence started impulsively from rest is 
considered, and will help clarify the issue of the Kutta condition. This is a classical 
problem (Prandtl and Tietjens, 1934). The airfoil is NACA 0012; the integral boundary- 
layer treatment is activated (see §8.2) and the Reynolds number is high. Figure 5 shows 
four snapshots of the flow. At first one sees a sheet of vortices leaving the leading edge and 
curling up. The starting vortex then grows in size. Since total circulation is conserved, an 
opposite circulation remains on the airfoil, contained in the boundary layers. A sheet of 
vortices connects the primary vortex and the trailing edge. This sheet clearly has a nonzero 


36 



circulation (velocity jump), since it undergoes a Kelvin- Helmholtz instability. Thus, the 
circulation on the airfoil is progressively brought up to its final value. In fact, the downwash 
of the vortex on the airfoil reduces the effective angle of attack, and tends to 0 only like 
1 It. The last snapshot is at much larger time; the wake now has zero circulation (note 
that it does not roll up any more) and the flow leaves the trailing edge smoothly. The 
force vector at early times shows a lower lift than in the final state, as well as some drag 
(induced drag from the starting vortex). At large times the drag is very small. 



Figure 5. Airfoil starting impulsively without initial circulation, 
a) tUoo/c = 0.02; b) tV^jc = 0.2; c) Woo/c = 0.3; d) tUoo/c — ► oo. 
c is the chord and the velocity. • • • vortices; ] force 


7 .6. Pressure and Force Extraction * 

The pressure was eliminated when the curl of (lb) was taken to derive (4) and is not 
necessary to advance the vorticity equation, but one often needs it as a diagnostic or for 
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the treatment of boundary layers. The wall pressure and the forces on the bodies can be 
obtained directly from the rate of creation of vorticity. Consider a wall with s and n the 
tangential and normal coordinates. Equation (7) can be written = — V.(Uu> — i/Vu>), 
i.e., the flux of vorticity is Uu> — i/Vui, which reduces to —uVu> at the wall since U = 0. We 
take the dot product with the normal vector n to obtain the flux through the boundary: 

— uduj/dn. This is the rate at which the wall is “creating vorticity”. Now consider the 
momentum equation (lb); at the wall it reduces to Vp = ^V 2 U or equivalently Vp = 

— v( V x u>) (using a vector identity). In 2-D, or in 3-D with a flat surface, one can show 
that n x (V x u>) + n.(Vw) = 0 (in 3-D, curvature seems to introduce another term). In 
those cases we then get n x Vp = udu>/dn , or in 2-D 


dp _ _ 0u> = a 2 r 
ds V dn ~ dsdi 


(39) 


The notation d 2 T/dsdt is appropriate since this is a rate of creation of circulation per unit 
length and unit time. This quantity is known, since we are creating the new vortices at 
every time step. It is then a simple matter to obtain the wall pressure in 2-D: equation (39) 
gives us the derivative of the pressure along the wall. We usually apply the trapezoidal 
rule and write 

_n±£ ki 

2At 


Pj+i ~Pj = 


(40) 


where the new vortices of circulation and T' + J are created during a step of length At. 
If we integrate (39) around the wall, since p is single- valued, we get 


</ — 

J dsdi 


= 0 , 


(41) 


that is, the total circulation created along each body is 0. This is how (38) was derived 
by Wu and Sankar (1980). Thus if we apply (38) the wall pressure given by (39) returns 
to exactly the same value after going around the body, which is not always the case with 
other methods for computing the pressure. A weakness of (39) is that it gives p only up 
to an additive constant. In many cases, the front part of the body is in contact with 
irrotational fluid from upstream; one can then obtain a good estimate of the additive 
constant by assigning the stagnation pressure to the front stagnation point (this neglects 
the unsteadiness in Bernoulli’s equation). Also, this method does not provide the pressure 
away from the wall. One could however use this wall pressure as boundary condition 
for the Poisson equation V 2 p = — V.(U.VU) (although the additive constants would be 
troublesome). 

The equation dpjds = d 2 T/dsdt which led to (40) is much less sensitive to the viscous 
details of the boundary layer (more specifically, the exact distribution of the vorticity in the 
normal direction) than dp/ds = —vduj/dn (used by Chorin, 1973). Thus (40) is very well 
adapted to vortex methods and an accurate wall pressure can be obtained even with a crude 
representation of the boundary layer itself. For instance if one has a steady attached flow 
the flux of vorticity along the wall is C/ 2 / 2 as mentioned earlier (§7.2). It is balanced by a 
flux of vorticity through the wall at a rate d(U 2 /2)/ds; thus we get dp/ds = — d(U 2 /2)/ds, 
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i.e., the well-known Bernoulli equation. Many procedures have been proposed to obtain 
the pressure from the time-dependent version of Bernoulli’s equation; these methods are 
very cumbersome in separated flows because of the numerous branch cuts in the velocity 
potential (one for each vortex). 

The global pressure force on a body is given by 


e 2 x 


(j> p(s) d\. 


After integration by parts and use of (39) this becomes 


e 2 x 



x ds, 


(42) 


(43) 


which is the negative of (14), applied to the newly created vorticity. This is not a coin- 
cidence. Neglect the viscosity for now. What happens is that the motion of the vortices 
under (26) preserves (14) (the momentum of the fluid), and that the creation of vorticity 
effects the exchange of momentum between body and fluid. Thus (43) amounts to a global 
momentum balance (but it can be applied to each body separately). One can check that 
(43) gives the right “apparent mass” for a body undergoing acceleration. A similar formula 
holds for the moment on the body, using the invariant (15) (Wu and Sankar, 1980). If 
the viscosity is not neglected, for instance if a random walk is used, then (14) is not (and 
should not be) conserved when a vortex collides with the wall, and this is how the friction 
force comes in. In practice, even without random walk there are such collisions due to 
noise in the velocity field near the boundary (see § 8.2 and 9.3 ) so that the pressure force 
and the momentum balance do not exactly agree, but this is a small effect, and (39) and 
(43) have been very satisfactory at high Reynolds numbers. 
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8. SEPARATE TREATMENT OF THE BOUNDARY LAYERS 


8.1. Overview 

The boundary layers are the most difficult regions for a vortex method for at least three 
reasons. First, the standard method uses radially symmetric elements and therefore is not 
adapted to the disparate scales near the wall (with finite-difference methods, high-aspect- 
ratio grids are routinely used near the wall). This results in significant noise in the normal 
component of velocity (recall figure 4b). Second, near the wall the viscous terms become 
dominant, and the available vortex methods do not treat them very efficiently. Third, the 
boundary-layer region is usually close to being steady in Eulerian coordinates: du >/dt is 
small, in contrast with the outer region where it is the Langrangian derivative Du>/Dt 
which is small. Thus Eulerian methods have an intrinsic advantage in the boundary layer. 

One response to the problem of disparate scales is Teng’s (1982) method, which uses 
elliptic vortices. Teng applied it only to a boundary layer, using identical ellipses and a 
random walk for the viscosity. He obtained good results, but apparently the method has 
not been used further. One could certainly extend the idea by using ellipses of different 
eccentricities; the vortices would be very flat near the wall and circular away from it. The 
ellipse eccentricity and orientation would be simple functions of the distance to the nearest 
wall. Note that such a concept is quite different from one in which the ellipse parameters 
would obey dynamical equations reflecting the straining and the self-induced rotation. The 
ellipses are introduced for kinematic reasons only; to obtain a more realistic flow field near 
the wall. A drawback of Teng’s method is that the formula for the velocity induced by 
an elliptic vortex is much more complex than for circular blobs, so that the method can 
be competitive only if an ellipse can replace at least several blobs. It would pay to find a 
simpler formula for an “elliptical blob.” 

Many methods can be found in the literature in which the no-slip condition is not 
enforced using the boundary elements as described in §7. Most often the no-penetration 
condition is enforced using images, and there are a few “separation points” at which 
vortices are created and released into the flow (Lewis, 1981; Stansby and Dixon, 1982). 
When the body has convex sharp corners, they are obviously going to cause separation, 
and it is consistent, to apply a Kutta condition. The rate of release of vorticity is U* / 2, 
with the appropriate sign, where U e the edge velocity of the upstream boundary layer. This 
prescribes the circulation of the new vortex, and the Kutta condition provides a condition 
on its position. Such methods satisfy the intuition and have the major advantage that 
only a few new vortices are created at every step, compared with hundreds in the method 
presented in §7.2. They are typical of applications on minicomputers. On the other hand, 
the use of the Kutta condition is still very delicate. Another weakness is that separation, in 
the sense of a transfer of vorticity from the wall region to the outer region, does not occur 
only at the obvious points. “Secondary separation” occurs at random on the back face 
of bluff bodies, and its neglect may be largely responsible for the excess circulation found 
by some studies and sometimes compensated for by arbitrary means (see §9.6). Stansby 
and Dixon (1982) emphasize the importance of secondary separation. In the long run, 
and especially for smooth bodies, methods based on a few selected separation points will 
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probably be superseded by those which allow separation in any region, depending on the 
boundary-layer dynamics. 

The idea of using an inviscid solver for the outer flow and a boundary-layer solver near 
the wall is old, and has been very successful for attached flows. In such a case one is 
primarily interested in the viscous friction drag, and possibly in the displacement effect 
on the outer flow. The trouble starts when the flow separates; the standard boundary- 
layer equations develop the Goldstein singularity. Any proposal to use the boundary-layer 
equations for a separated flow should include a discussion of this singularity. It is all too 
easy to mistake its effects for purely numerical problems, or conversely to suppress them 
using numerical dissipation. 


8.2. Simple Delay of Separation 

As mentioned already the vortices induce a fair amount of noise near the wall, which 
tends to scatter them. As a result, some move away from the wall, and others collide with 
it. The ones that collide can be pushed back, or they can be absorbed and their circulation 
be reintroduced a distance <f 0 from the wall (see §P.5). Statistically, the vorticity tends to 
drift away from the wall. As soon as there is an adverse pressure gradient, the straining 
adds to this effect, and separation occurs. This is a mild case of the situation described for 
sharp edges in §7.^. The pure vortex method consistently predicts separation at the onset 
of an adverse pressure gradient, which is essentially the behavior of a laminar boundary 
layer. With turbulent boundary layers, separation can take place much farther along the 
wall and a 2-D method is unable to reproduce this by itself. 

This led to the idea of predicting the separation point by a separate method, which can 
take the Reynolds number and the effects of turbulence into account, and intervening in 
the dynamics of the vortices to delay their separation up to the correct position. If this 
intervention is performed in a sensible manner, we can retain the smooth blending between 
the attached boundary-layer vorticity and the outside vorticity (this is not achieved with 
the methods that release vortices only at a few points and do not enforce the no-slip 
condition). 

At every time step, the pressure distribution is computed using (40). With it, the two 
boundary layers starting at the front stagnation point of each body are calculated. So far, 
only integral methods (Thwaites’ and Head’s, see Cebeci and Bradshaw, 1977) have been 
used, in conjunction with Granville’s transition criterion, and the time-dependent effects 
on the boundary layers have been neglected. Thus, there is room for experimentation 
with more sophisticated boundary-layer solvers. The only information needed from the 
boundary-layer solution are the separation points (we do not attempt, to continue the 
boundary-layer solution beyond separation). On each body, downstream of these two 
primary separation points the vortices are free to behave normally, so that secondary 
separation can take place. 

Upstream of the separation points indicated by the boundary-layer solution, we prevent 
the statistical drift and separation of the vortices simply by re-absorbing the new vortices 
after only one step. A fresh regular row of vortices carries the boundary-layer vorticity 
as in figure 4, representing a thin attached boundary layer. Just beyond the indicated 
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separation point, the vortices are suddenly free and subjected to the straining associated 
with the adverse pressure gradient, and they leave the vicinity of the wall. This procedure 
achieves the goal of making the vortices separate approximately in the right place. 



* 



Figure 6. Simulation of the flow over an airfoil, a) without delay of separation; 
b) with delay. ••• vortices; f lift and drag 


Figure 6 shows the flow over an airfoil calculated without and with separation delay. In 
figure 6a, observe how the vortices remain attached on the lower surface, thanks to the 
favorable pressure gradient, but separate just downstream of the leading edge on the upper 
surface, in a mild adverse gradient. This is what laminar boundary layers would do. In 
figure 6b with a very high Reynolds number the integral boundary-layer solvers predicted 
turbulence, and separation only at the trailing edge. This activated the separation-delay 
device and we are left with attached boundary layers and a thin, smooth wake. The lift is 
higher, and there is almost no drag. Separation delay was also applied in figure 1; it made 
a difference on the upper surface of the slat, and on element 4. It was fully responsible for 
the Reynolds- number dependence that was described (§ ^ . J). 

This device has been used extensively, and is very helpful in many cases (figures 1, 5, 
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6b, 7, 13 and 14), but it also has weaknesses; it cannot account for reattachment. The 
weakest part of the solution is the transition prediction (but this is true of any method). 
To obtain a flow pattern that agreed with experiment in figure 7, the transition criterion 
had to be adjusted; the increase in R 6 from the instability point to the transition point was 
taken as half of Granville’s correlation. Another minor problem occurs when one cannot 
find the front stagnation point because the flow is too disturbed in front of the body. This 
occurs in the rotating-stall calculation, see figure 14. In that case, one just bypasses the 
boundary-layer solution for this body at this step. Finally, the integral methods have a 
tendency to predict separation a few wall points downstream of sharp corners, and one 
needs to override this prediction and set the vortices free a few points upstream of the 
corner. Because of these issues, when starting with a new shape it is good practice to first 
simulate the flow without separation delay, then activate the device if there is a need for 
it, and carefully check the behavior of the boundary layers. 

Figure 7 shows a case in which separation delay made a significant difference. An NACA 
0012 airfoil is oscillating from 5° to 25° at a reduced frequency of 0.25 (Spalart and Leonard, 
1981). The overall flow strongly depends upon separation on the upper surface just behind 
the leading edge, and this separation is controlled by transition in the boundary layer. The 
Reynolds number was 2.5 10 6 . The forces and moments on the airfoil were in very good 
the agreement with experiments (Spalart, Leonard, and Baganoff, 1983). 

8.3. Chorin’s Tile Method* 

In 1978 Chorin introduced a vortex-like method for the boundary layer. As before, the 
elements carry circulation, move with the fluid, and undergo a random walk for viscous 
diffusion. However they are now “sheets” or “tiles” and the velocity is computed in a much 
simpler manner thanks to the boundary-layer approximations (the tiles interact only with 
their neighbors). This method is much cheaper computationally than the usual method 
or Teng’s method (1982), and is designed for use in a hybrid algorithm, with conventional 
vortices away from the wall. Circulation is exchanged by the two subdomains, primarily 
in the direction sheet — ► blob. Chorin discusses the Goldstein singularity and the need 
for a viscous-inviscid coupling strategy to remove it, but the subsequent users of the tile 
method completely overlooked this issue. 

Cheer (1983) applied the complete hybrid method to the flow past a circle and Joukovsky 
airfoils, with Reynolds numbers of about 1000. She obtained flows with “laminar” behavior, 
i.e., separation as soon as the boundary layer enters the adverse pressure gradient, and 
satisfactory drag and lift values. Ghoniem and Gagnon (1987) applied it to an internal 
flow at rather low Reynolds numbers, which seems a little superfluous since there were no 
thin boundary layers (the inflow was a parabolic profile); quite probably a pure vortex- 
blob method would do just as well. In summary, the tile method succeeded in introducing 
the boundary-layer character for the near-wall flow; it is fast and solves the first problem 
described at the beginning of §8 (the scale disparity), but not really the other two. The 
convergence of the random-walk method is rather slow here as it is with blobs. Also, the 
Goldstein singularity is a major problem and Chorin himself (1978) did not assert that it 
was fully solved. 
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8-4 ■ Coupling with a Finite -Difference Solver * 

As mentioned above, Eulerian methods have significant advantages for boundary layers, 
especially if they are steady or quasi-steady, which makes the time integration undemand- 
ing. Furthermore, the location of the steep gradients is known and it is easy to cluster the 
grid points near the wall. This motivated an effort to couple an inviscid vortex solver, for 
the outer region, and a finite-difference solver, for the wall region (Spalart, Leonard, and 
Baganoff, 1983). 



Figure 8. Sketch of the hybrid finite-difference/vortex algorithm, o vortices; — 
grid lines; — > velocity 


The inner region is a thin strip along the wall, see figure 8. The boundary-layer vor- 
ticity equations are solved using a standard, implicit, centered finite-difference scheme (T. 
Pulliam, personal communication, 1981). Simple transition and turbulence models are 
used. The challenge is in the coupling with the outer region. The two regions exchange 
circulation as in the tile method, according to the direction of the velocity at the interface. 
The coupling of the velocity fields is more delicate. A viscous-inviscid coupling was imple- 
mented based on the displacement thickness (the centroid of the inner-region vorticity). 
The procedure is intuitively correct and very similar to the ones used in attached flows, 
but there is no proof that it actually removes the Goldstein singularity, especially when 
large amounts of vorticity cross the boundary. Some numerical difficulties (oscillations) in 
the separation region suggest that it does not (see Spalart et al., 1983, for details). Partly 
for this reason, this method has not been as robust, and has received much less use than 
the simpler “delay” method (§ 8.2 ). 

The method however gave encouraging results for the flow past a circle. The “drag 
crisis” in this flow at Reynolds numbers near 3 x 10 5 is one of the major challenges in 
the field of viscous flows. The flow was calculated for Reynolds numbers between 10 4 
(the lower limit for thin boundary layers to form) to 3 x 10 7 (beyond which there were 
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Figure 9. Drag coefficient of a circular cylinder. — average of experiments; 
o calculation 

no experiments). Figure 9 shows that a decrease in drag coefficient was observed, but 
that the crisis itself (including the sudden drop to very low drag values and the steady 
lifting states observed in experiments) was not. The drag decreased at higher Reynolds 
numbers because the boundary layers transitioned; the turbulence model then kept the 
vorticity attached longer, so that it crossed the boundary into the outer region farther 
downstream. This is illustrated in figure 10; the mean streamlines, shown at Re = 10 s and 
10 6 , reflect the position of the separation region. By symmetry, only half of the flow needs 
to be shown. The Reynolds-number effect was qualitatively correct, and the quantitative 
agreement with experiments, over a wide range of Reynolds numbers, was fair. 

For the future, the use of an Eulerian solver for the inner region should be the most 
efficient procedure. On the other hand, it is not the most elegant since it is a hybrid, 
and the coupling and singularity issues definitely need to be re-examined. One could of 
course solve the full Navier- Stokes equations in the inner region; the obstacle there is 
that one would need to solve a Poisson equation to obtain the velocity field, instead of just 
integrating the equations u y = — and v y = —u x up from the wall. The velocity boundary 
condition on the interface would be a serious problem. 
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9. PRACTICAL ASPECTS 


9.1. Time-Integration Schemes 

A wide variety of time-integration schemes are used in computational fluid dynamics, 
and the choice is often a matter of taste in addition to the technical concerns. Still, one 
can narrow down the choice considerably in the case of vortex methods. The first point is 
that all the schemes used have been explicit schemes. To use an implicit scheme one would 
need to invert (26) or at least a linearized form of it, and (26) is too complex for this to be 
done efficiently. Moreover the usual stability advantage of implicit methods is much less 
crucial with vortex methods than with grid-based methods (see below). The second point 
is that the most costly part of the algorithm, by far, is the evaluation of the time derivative 
using (26). Also, the calculations use a relatively small amount of memory, so that storing 
old values of the time derivative is not a problem. This suggests that multi-step schemes 
are preferable to predictor-corrector schemes. 

Based on these remarks the author has used the Adams-Bashforth second-order scheme, 

xj{t + At) *Xj(t) + A/Qu(x j ,f) - |u(x,-,< - A t)^j (44) 

for the integration of (20), which represents a good compromise between accuracy, sim- 
plicity and cost. As usual, the first step for each new vortex is done using the Euler 
scheme. A majority of the studies in the literature used either the Euler first order scheme 
or Runge-Kutta schemes of second or fourth order (the Runge-Kutta second-order scheme 
is also known as “modified Euler” or “Huen’s scheme”). The use of the low-order Euler 
scheme, when the Adams-Bashforth scheme is hardly more costly, was motivated either by 
a question of consistency with the random-walk algorithm (in Chorin’s method) or by a 
faulty concept of numerical dissipation, as discussed below. 

The integration of (20) is most difficult when large accelerations are experienced by a 
vortex, typically when it is close to another vortex. Therefore a good model problem is the 
motion of just two vortices (Spalart and Leonard, 1981). The problem can be expressed 
in complex notation and simplified to 


dz ift 
dt z * ’ 
2 ( 0 ) = 1 


(45a) 

(456) 


where z is the separation vector between the two vortices, z* is its conjugate and ft is a 
frequency. The exact solution is z = exp (iftf): the vortices orbit around each other with 
angular velocity ft. Their distance is constant. It is a good exercise to program (45) and 
observe the behavior with different numerical schemes. 

One finds that with any scheme, if the time step At is large compared with 1/ft the 
solution is inaccurate. However due to its nonlinearity (45a) reacts very differently from 
the linear equations that are usually studied, e.g., the equation dz/dt = iftz which has 
the same exact solution. Usually, as an explicit scheme becomes inaccurate (because At 
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is too large) it also becomes unstable and the solution grows to infinity (it “blows up”). 
With (45a) the distance between the vortices grows until |z| 2 /fl is of the order of At, 
then settles down. The vortices now orbit at a lower angular velocity, fl/|z| 2 . With some 
schemes |z| will shrink again and in some cases one gets an oscillation, especially with 
multistep schemes as the primary and spurious roots alternatively dominate (Spalart and 
Leonard, 1981). With the Adams-Bashforth scheme, such unphysical behavior has not 
been observed. 

The lesson to be learned here, which holds in practice with many vortices and has 
been recognized by many authors (Milinazzo and Saffman, 1977, etc.), is that in almost 
all situations time-integration errors tend to scatter the vortices. The scattering occurs 
whenever the local rotation rate fl is large compared with 1 /At and corrects itself, since 
the scattering precisely reduces fi. The corollary is that a sustained instability, with the 
solution “blowing up” to infinity, never occurs. In that sense, the method is “robust”. 
However, a solution can remain bounded and be very inaccurate. Also, a scattering of 
vorticity reduces the velocity fluctuations and the kinetic energy. In that sense, it represents 
a numerical dissipation. Lagrangian methods, unlike fixed-grid Eulerian methods, allow a 
perfect convection of flow structures by a uniform velocity field; however when gradients 
are present inaccuracies show up as with any other method. 

The scattering of vortices has implications for the invariants (11, 12, 15, 16). The total 
vorticity (9, 13) and the impulse (10, 14) are linear quantities in terms of the vortex 
positions Xj. Therefore with the common time-integration schemes which are all linear, 
if (9, 10, 13, 14) are conserved by the spatial discretization they are still conserved even 
allowing for time-integration errors in (20). This is not so for the nonlinear invariants 
(11, 12, 15, 16); they are only “semi-conserved” and have often been used as measures 
of the integration errors. For example Nakamura, Leonard, and Spalart (1982) defined 
a global measure of the numerical dissipation, an “effective viscosity,” based the time 
evolution of (15), primarily to verify that it was purely caused by time-integration errors 
(indeed, it was). Some researchers pushed this reasoning one step too far, concluding that 
time-integration errors were desirable, since they introduce the viscosity that is so difficult 
to introduce by other means! They then use the most inaccurate and dissipative (when 
applied to (45)) scheme, the first-order Euler scheme. Nakamura ct al. never implied that 
the effective viscosity had any meaning locally. Moreover, if errors are what one is after, 
it is cheaper to just increase the time step, rather than switching to a lower-order scheme. 

One final remark about steady flows. If the velocity field is naturally time independent, 
grid-based methods can usually find the solution faster by direct solution or by relaxation 
techniques than by a straightforward integration in time. Vortex methods do not have 
this option and treat any flow as time dependent; the vortices move even in a steady flow. 
Fortunately, most separated flows are strongly unsteady, so this difference is not of much 
importance in practice. 


9.2. Core Functions 

In a separated-flow calculation one will need to compute F and F many times, so that 
it is desirable for both to have simple formulas. Consequently it is convenient to choose F 
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first. As mentioned earlier the Gaussian core (21) (but with a fixed radius a substituted 
for \/ui) is attractive but is expensive to compute, and yields only second-order accuracy. 
A very simple core profile was used in all the cases shown in these notes: 



(46) 


It is second-order accurate too. Its only weakness is that the vorticity decays rather 
slowly at large distances, like r~ 4 , so that the residual vorticity mentioned in §7 .3. is not 
negligible. A vorticity distribution with bounded support or at least faster than r~ 4 decay 
would be preferable as far as the interaction with walls is concerned. Hald (1979) and 
Beale and Majda (1982, 1985) give formulas for higher-order cores. Recall that Perlman 
(1985) found that these performed better only if the vortices remained well-ordered, so 
that in complex separated flows their advantage is probably marginal. 


9.3. Control of the Vortex Count 

In most applications, many new vortices are created at every time step and added to 
the list. For instance in the calculation of figure 1, there are about 1300 vortices in total, 
and 249 new ones are created along the wall at each step. Clearly, useful calculations are 
possible only if vortices are deleted at about the same rate they are created. 

First take the simple case of the mixing layer of figure 2. At each step, one new vortex 
is created at the trailing edge, and one “old” vortex is deleted near the outflow boundary. 
The easiest procedure is to delete the oldest vortex; a slightly better one is to delete the 
vortex that is farthest downstream. This makes little difference, because the last 10 or 20% 
of the domain are not included in the “measurements” (e.g., velocities, Reynolds stresses) 
anyway since that region is disturbed by the lack of vortices beyond the outflow boundary. 
Either method allows one to generate a satisfactory statistically steady state of the flow, 
and to maintain it as long as desired. 

In the less simple cases with solid bodies, two operations allow one to control the number 
of vortices. The first one is the absorption of vortices by the wall. When a vortex collides 
with the wall, one has several options. One can reflect it as in an elastic collision, or 
just push it back to the wall. A third option which is consistent with the treatment of 
viscous boundaries described in § 1.2 is simply to delete the vortex. If this is done before 
the evaluation of t/’ y for (31), the right-hand side of (31) will reflect the loss of circulation 
near the wall, and the new vortices will make up for it (also adjust the last line of the 
right-hand side, (32)). This absorption procedure has the advantage of helping to contain 
the number of vortices, and of slightly reducing the noise near the wall. 

The absorption of vortices is helpful, but is not enough since it is not under the user’s 
control; a majority of vortices still leave the vicinity of the wall and move into the wake. 
Thus one needs a more active and adjustable means of limiting N . This can be done by 
merging vortices. The method used is adapted from Rogallo’s (personal communication, 
NASA Ames Research Center, 1979). Deffenbaugh and Marshall (1976) also used a simple 
merging device. At every step, pairs of vortices are merged into single vortices. The goal 
is to achieve this with the least amount of modification of the velocity field. 
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Since a large local modification is unavoidable, it makes sense to require that the velocity 
in the far field be altered as little as possible. Suppose the vortices j and k are merged, 
and the new vortex is given the circulation T* and the position z\ Using complex notation 
and omitting some constant factors, the change in the velocity at z is: 

r Tj Tj, 

z - z' z - Zj z - z k 


(47) 


The Taylor expansion for large \z\ is 


r 
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(48) 


The two leading terms are removed by taking T' = Tj + T k and z 1 = (TjZj + T k z k )/T': 
the new vortex inherits the sum of the circulations of the old ones, and is placed at their 
centroid. This could be expected; the circulation (13) and the vortex centroid (14) are 
conserved. On the other hand, (15) and (16) cannot be conserved (recall that they were 
already not conserved because of time-integration errors). 

The next term in (48) cannot be removed and we take it as the merging criterion, i.e., 
we merge a pair only if this term is below some tolerance. After insertion of T ( and z' and 
taking of the modulus it becomes: 


ir,r fe [ \ Z j - z fc j 2 
|r; + r fc | |*|* 


(49) 


Consider the two factors in this estimate. The first one encourages merging of weak vortices 
(small Tj and/or IT), which is sensible since weak vortices cost as much as strong ones, 
but carry less information. It also discourages merging of vortices with nearly opposite 
circulations (small |Tj + T*.|) which also makes sense since in such a merging z' ends up 
far away from Zj and z k which is not natural. 

The second factor encourages merging of close vortices (small \zj — z k \) as one would 
expect. We still need to prescribe |z| in the denominator. The most sensitive place in the 
flow is the wall region; therefore we insert the distances dj and d k from the vortices to the 
closest wall. The formula that was finally adopted is 

!££*] h - *>1* 

|r, + r*| (B„ + + i k yi' 

Notice how the formula was symmetrized with respect to dj and d k . The length parameter 
D 0 allow one to shift resolution closer to the walls (small Do) or into the wake (larger D 0 ), 
but the calculations show little sensitivity to Do over a wide range (Spalart, Leonard, and 
Baganoff, 1983). A typical value is 10% of the chord. 

Finally, the quantity V 0 in (50) has the dimensions of a velocity and is the tolerance. 
Typical values are small: less than 10 -4 times the freestream velocity. In practice it is most 
convenient to let the program adjust V 0 during the run. One chooses a target number of 



52 


vortices, and employs a feedback mechanism that raises V 0 if there are too many vortices, 
and vice-versa. The merging procedure described here has been very convenient in practice, 
is intuitively correct and has a rational basis in the Taylor expansions. On the other hand, 
merging is a very noisy operation locally, and its effect on the convergence is unknown. 
Therefore, it is not appropriate for all applications. 


9.f Efficient Programming; Operation Counts 

Vortex codes are relatively easy to optimize, because almost all the operations are spent 
computing the A r 2 interactions in (26). Thus one really needs to optimize only about ten 
lines of code. There are a number of nearly-obvious tips. The division by 27 r and the x 
product should be done outside the 0(N 2 ) loop. The formulas for the two interactions, 
k — ► j and j — > k, have a lot in common, namely the computation of 

( x j ~ x fe) r / l x j ~ x fe 

|Xj-Xfc| 2 \ <7 

After this, the circulations T & and Tj are presumably different and the overlap ends. One 
should have an outer loop ( k = 2 ,N) and an inner loop (j = l,k — 1) and compute (51) 
only once. There is no need to store it (which would occupy N 2 words of memory); use it 
right away to increment Ufc and Uj (or store all the ji — > fc terms and sum them up outside 
the j loop but inside the k loop). In 2-D there is no need to compute any square root, 
and one can avoid exponentials. The operation count is 8 N 2 with the simplest core, (46), 
and would be hardly higher for more complex but algebraic cores. In 3-D, the symmetry 
between k — > j and j —>■ k is useful too. One should also precompute quantities like the 
tangent vector outside the 0(N 2 ) loop. 

In practice, periodic calculations could be much more expensive than unbounded ones, 
because the stream function and its derivatives involve the computation of several com- 
plex exponentials N 2 times (see §10.1), and exponentials cost many floating-point op- 
erations. However, this drawback can be completely eliminated with a little extra pro- 
gramming. One needs the values of exp(iTr(zk — Zj)/p) for all values of k and j. What 
one does before entering the 0(N 2 ) loop of the program is to precompute and store 
exp(i7TZfc/p) and exp( -inzk/p) for every k. Then the operation that is done N 2 times 
is exp(iirzk/p) x exp(—iirzj/p): the exponentials have been eliminated. This can easily 
speed up the execution by a factor of 4 or more, depending on the computer. In fact, 
the 0(N 2 ) component of the operation count is only 18iV 2 so a well-programmed periodic 
method is only about twice as expensive as the unbounded-flow method. Thus there is no 
reason to approximate (52) by taking the first few terms in the series, as has sometimes 
been done in the past. The same trick applies to the computation of the stream function 
at the wall points, for (31). 

Another tip concerns the use of complex variables, in 2-D. They are convenient, but 
degrade the efficiency of some computers to some extent (the machines that slow down 
when the loops have skips of 2). Some machines even refuse to vectorize any loops involving 
complex variables. This should be checked carefully before writing a code. 
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9.5. Flow Charts 

Here we give flow charts of two programs, named KPD1 and KPD2. Both are available 
from COSMIC (a part of NASA) in the U. S., as are their periodic counterparts for 
cascades. A listing of KPD1 also is included in the author’s thesis (Spalart, Leonard, and 
Baganoff, 1983) (there is an error on line 15, p. 107; replace NW by NWALL(L)). KPD1 
is the basic bluff-body code with no boundary layer treatment; KPD2 has the integral- 
method treatment of the boundary layers (the “separation delay” of §#.2). The flow charts 
provide an overview of the material covered until now, and specify the order in which the 
various operations are performed. 



Figure 11. Flow chart of the KPD1 code 


9.6. Assessment of the Accuracy 

This can be delicate in practice. Some numerical methods (in particular with explicit 
time-integration schemes) give the user a strong warning when the time step is too long: 
they “blow up”. Other methods (spectral, for instance) generate “wiggles” which are also 
a warning that finer resolution is needed. Vortex methods do not give such warnings; in 
most applications the vortices look disordered, but the calculation never blows up. This 
is partly due to its automatic conservation of some key integral quantities (9, 10, 13, 14). 
Again, the possibility that the calculation is not accurate on small scales but is accurate 
on the intermediate and large scales is likely. 

The most prominent result of a bluff-body calculation is the average drag coefficient 
Cd • We treated four shapes with sharp corners, so that the primary separation points 
were fixed, and compared with experimental results (Spalart, 1984, and unpublished). 
The shapes were a square, an equilateral triangle in two positions and a vertical flat plate 
or thin airfoil. They were all treated with the same code and roughly the same resolution 
(typically 1000 vortices). The results are summarized in a table (the — * indicates the 
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Figure 12. Flow chart of the KPD2 code 


direction of the freestream velocity). 


Shape 

Cx?(calc.) 

Ci)(exp.) 

— > □ 

1.5 

2. 

— > > 

2.9 

2. 

— > <1 

2.4 

1.3 

- 1 

2.8 

1.75 


The uncertainty is at least ±0.2 on all the values. Even then, the disagreements are 
significant. Note also that the calculation often over-predicts Co, hut not always. In 
general, the shedding frequency (Strouhal number) is better predicted than the drag. These 
results are not very satisfactory, but we were unable to improve them or to convincingly 
demonstrate the convergence of the method, even via large variations of the numerical 
parameters (N , At, d 0 , etc.). The above table provides a severe but conservative estimate 
of the performance of the author’s method. In some cases, the agreement with experiment 
is much better, within 5% (see §10.1). 
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Many of the comparable methods also overpredict the drag. One reaction to this is to 
introduce a “fudge factor” to bring the drag down. Arguing that cancellation of vorticity 
of opposite signs occurs in reality but not in the calculation (which is not even obvious) 
some authors arbitrarily make the circulation of the vortices decay as they move into the 
wake. By adjusting the decay rate one can get any drag one wants. In the author’s mind 
the idea of invoking Kelvin’s theorem, building a method around it, and then violating it 
just to fine-tune one global quantity, C'c, is strange. 

Milinazzo and Saffman (1977) already pointed out that “comparison with gross experi- 
mental features, such as drag, is not conclusive, and in any case the real flows are turbulent, 
and therefore three-dimensional, at large Reynolds numbers.” This point is valid, espe- 
cially if the experiment al flow has large-scale three-dimensionality. It is possible that much 
of the drag difference is due to this effect. The author’s calculation of dynamic stall (figure 
7) agreed with experiment better than those for steady bodies. In this case, the whole 
span of the wing oscillates in phase, which makes the flow more two-dimensional. The 
degree of three-dimensionality of the flow past 2-D shapes, such as a circular cylinder or 
a backward-facing step, is a matter of debate and the subject of ongoing experiments. 
Reliable numerical simulations should also be obtained in the next few years. 

In summary, there is no recipe for the estimation of the errors in a vortex calculation of 
a complex flow like figure 1. It is essential to monitor the vortices and velocity field as in 
figure 1 or 4b, but this is only as good as the user’s intuition of what the flow should look 
like. Finally, it should be remembered that we are looking at very challenging flows. If one 
tackled the multi-element airfoil with an existing finite-difference method, very probably 
the best one could do is a rather crude calculation which would not enter its convergence 
pattern, at least not with a reasonable number of points. One would also be dependent on 
crude transition and turbulence models similar to those we used, so that the viscous effects 
are very delicate with any method. When applied to complex flows, the vortex method 
will not yield perfect results, but it will yield very useful ones. 
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10. TREATMENT OF CASCADES AND SCREENS 


10.1. Spatially -Periodic Flows 

The computation of 2-D flows that are spatially periodic in one direction is almost as 
simple as that of unbounded flows, and is useful for the treatment of temporally-developing 
mixing layers and of cascades, for instance. Complex notation is convenient here. Lamb 
(1945, p. 224) gives the formula for the stream function generated by a row of vortices of 
circulation T at the points zo + np, n being any integer ( zo and p are complex): 



7r(z - Z 0 ) 

P 



(52) 


Compare this formula with (8). This function can be regularized by a “core” as in un- 
bounded flows and a convenient one is the analog of (46): 
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The rest of the method is identical to the unbounded case; one simply substitutes (53), for 
the stream function, and its derivatives, for the velocity. See for a tip in programming 
(53). To make the visualizations consistent one chooses a “main period” between some 
value z\ and zj + p; if a vortex moves out of the main period (i.e., if the real part of 
(z - z\)/p leaves the interval [0,1]) the vortex is translated by ±p so that all the vortices 
can be seen next to each other. 

An example of a periodic calculation is the study of a row of chevrons, shown in figure 13. 
These chevrons were candidates for an acoustic barrier around the inlet of a wind tunnel at 
NASA Ames Research Center. The pressure drop across the barrier was estimated to be of 
the order of 5 to 10 times the dynamic pressure, based on the ratio of the open area to the 
pitch of the cascade (essentially the argument of Batchelor, 1967, p. 375). The calculation 
predicted a much larger value: 86 times the dynamic pressure. The issue was settled 
when experiments by Professor J. Foss (Michigan State University) produced a value of 
82. Such losses were unacceptable, and the design was dropped. The open-area argument 
vastly underestimated the loss of energy in the flow due to the repeated separations. This 
was one of the more successful and directly useful applications of the method. 


10.2. A Simple Model for Screens* 

Screens are often used in wind tunnels to generate homogeneous turbulence, to introduce 
a pressure drop (e.g., for mixing-layer studies), or sometimes to smooth out an irregular 
flow and prevent separation on a cascade. Here we present a simple, macroscopic model 
that can account for the latter three effects and is a good exercise in itself. It was used in 
the chevron calculations (see figure 13). It could be used for nonperiodic geometries too. 

Consider for simplicity a screen along the y axis in 2-D. We model it as a pressure jump 
A p(y) = p(0 + , y) — p(0~ ,y). The effect of the screen on the fluid is equivalent to a body 
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Figure 13. Flow past two rows of chevrons. • • • vortices; — » forces; — streamlines; 
screen 


force f = A p(y)f>(x)e x with 6 the 1-D Dirac distribution and e x the unit vector in x. To 
obtain the effect on the vorticity we take the curl of f which is ~d(Ap)/dy 6(x). There is 
a generation of vorticity at the rate —d(Ap)/dy per unit length in y and unit time. This 
is easily modeled by creating new vortices along the screen at each time step. 

We still have to prescribe the value of A p. In the chevron calculations a pressure drop 
proportional to the local dynamic pressure was assumed, that is: A p = — C u\u\/2 where u 
is the local velocity in the x direction and C is a positive constant and in general depends 
on the solidity of the screen, the shape of the bars, and so on. C was set to 2. Near 
the center of figure 13 there is a jet (closely spaced streamlines, corresponding to high 
velocity); notice how it spreads out when it encounters the resistance of the screen, as 
expected, resulting in a slightly more uniform flow behind the screen. 

So far, the simplest method has been used for the time integration: the value of u at 
one step is used to compute A p and the circulation of the new vortices for the next step. 
This method worked with the parameters used in the chevron study, but will produce an 
instability for large enough values of the product C\u\At/ Ay, where Ay is the spacing of 
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Figure 14. Propagation of a stall cell in a cascade. • • • vortices; f forces; 
— streamlines 


the points that describe the screen. For large values of the product (in the limit C — * oo 
one actually has a solid wall) one will need a fully coupled algorithm, instead of one which 
goes back and forth between the fluid and the screen. 


10.3. Computation of Rotating Stall in a Cascade* 

The periodic version of the method was also applied to 2-D cascades of airfoils with 
the hope of observing “rotating stall”, or more appropriately “traveling stall” (since the 
cascade is actually not. circular). When a compressor stalls, the stall pattern is often 
uneven from blade to blade, and stall cells tend to propagate around the compressor at a 
fraction of the velocity of the blades. Such a flow involves a multiply-connected domain 
(since several blades are involved) and large regions of separation, which makes it a good 
candidate for the vortex method. 

Figure 14 is taken from Spalart (1984, 1985). Five NACA 0012 airfoils are included 
within the numerical period p. They are at 45° to the plane of the cascade. The incoming 
flow, of magnitude Uo, gives them a 30° incidence. Let c denote the chord of the airfoils; 
the pitch was chosen equal to c. The time interval between frames in the figure is 2c/Uq- 
The computation was started from irrotational flow with a dipole in the wake (to break the 
blade-to- blade periodicity) and has been running long enough to reach a well-developed 
state (see Spalart, 1985, for the transient behavior and other values for the stagger and 
inflow angles). One observes a stall cell which alters the flow near two of the blades, 
and drastically reduces or even reverses the mass flux between them. The cell propagates 
upwards at about 40% of Uo- This figure, and the flow pattern, are in qualitative agreement 
with flow visualizations, but more detailed comparisons were not attempted. The method is 
limited by its reliance on rather crude transition and turbulence models (as in the dynamic- 
stall study), by the resolution which is becoming marginal (100 wall points and 300 vortices 
per blade), and by the fact that it is 2-D. However Speziale, Sisto, and Jonnavithula (1986) 
extended the study to cambered airfoils and a wider parameter range. 


11. TIME-SAVING ALGORITHMS 


11.1. Vortex-in- Cell Methods* 

These methods discard the Green ’s-function approach (3, 6, 8) and solve the Poisson 
equation for the stream function by a grid method (Christiansen, 1973). Typically, with 
JV* grid points the solution is obtained in 0(N' log(JV')) operations, so if N' is comparable 
with N the vortex-in-cell method is much faster than the 0{N 2 ) method. In practice N' 
is at least as large as N , and could be much larger if the vorticity distribution is very 
intermittent. The procedure is to define a grid and to obtain vorticity values at the grid 
points according to the circulation of the vortices (if any) present in the adjoining cells. The 
vorticity is “transferred” to the grid, but only for the Poisson solution; the vortices retain 
their identity. The beauty of this method is that it still has the low-dispersion advantage of 
other vortex methods. Once the stream function is available it is differentiated to obtain 
the grid values of the velocity. The velocity is then interpolated to the position of the 
vortices and the vortices are moved, thus completing the work for this time step. 

The best applications are those in which N' can be kept to a value close to N , that is, 
when the vortical region is rather predictable and of simple shape, for instance an internal 
flow or a mixing layer (Aref and Siggia, 1980, Couet, Buneman, and Leonard, 1981). In 
such cases, one can save a large fraction of the computational effort. The method is not 
grid-free, and since boundary conditions are needed on if> the grid has to be body-fitted. 
This would be a serious drawback in engineering applications. Note also that the transfer 
and interpolation steps have a cost only O(N), but. are hardly vectorizable. Another topic 
of research is in the local errors associated with the exchange of information between the 
regularly-arranged grid points and the vortices, which can fall in any region of the grid 
cell. The grid destroys the isotropy of the algorithm, and in general a vortex will have a 
spurious self-induced velocity. See the work by Couet. et al. (1981) and Anderson (1985b). 


11.2. Lumping Methods * 

These methods still use the Biot.-Savart approach, but reduce the operation count by 
approximating some of the interaction terms in (26), using Taylor expansions. Spalart and 
Leonard (1981) presented such an algorithm, claiming an operation count, of order JV 3 / 2 . 
L. Greengard and Rokhlin (1987) claim O(N) for their algorithm. Divide the domain 
where the vortices are into a regular array of cells. The cells are of size /, large enough to 
contain many vortices each (in contrast with the vortex-in-cell method). They can overlap 
with the solid body; in that, sense the method is “practically” grid-free since the grid is 
not body-fitted. Now consider two cells that are separated by at least, a few times l. The 
interactions are proportional to 1 /z, using complex notation. This function is analytic and 
has a rapidly converging Taylor expansion. If Zl and are the cell centers, Z{ is in cell 
L, and Zj in cell K , we Taylor-expand the function 1/z in the vicinity of Zk — Zl; see 
figure 15. 

There are several options. The requirement is to incur a relative error of less than some 
small number e on any of the terms. Spalart and Leonard chose to use a varying number 


PRECEDING PAGE BLANK NOT FILMED 


61 




Mimmmx &ahk 




Figure 15. Schematic of cell-to-cell interactions in the “lumping” methods 


of terms depending on the value of \Zjc — Zl\/1 (more terms if it is small). If the cells are 
neighbors the interactions are computed vortex-by-vortex, because the Taylor expansion 
fails (i.e., it is too close to its radius of convergence). One can optimize the number of 
cells, and the best value is of order N 1 ? 2 . The overall cost, for a fixed value of e, is then 
of order N 3 / 2 . Greengard and Rokhlin (1987) carry the idea further by nesting the cells. 
There are many levels of lumping. In essence they keep | Zjc — Zl\/1 roughly constant, 
instead of varying the number of terms. Spalart and Leonard (1981) wasted terms on the 
very- well separated cells (large \Zk — I/O- The nesting removes the need for computing 
many vortex-by-vortex terms. On the other hand, one needs many more cells. Greengard 
and Rokhlin’s arrive at an operation count of order N log 2 (e). 

Greengard and Rokhlin’s (1987) statement that the cost is of order N is not rigorously 
correct, because e is not a constant (there is no reason why it should equal the precision 
of the computer). For the method to converge, c must tend to zero like some power of 
N as N tends to oo. Thus the cost is really of order N \o^{N). Similarly, Spalart and 
Leonard’s N 3 / 2 estimate should be revised upwards, at least to IV 3 / 2 log(TV), because as 
e — > 0 one needs to include more Taylor terms for a given separation. A factor log(JV) 
is not very important especially since larger values of N improve the throughput of the 
computer, but one likes to be rigorous if possible. 

In practice, Spalart and Leonard’s (1981) algorithm was helpful on serial computers 
(CDC 7600) with e = 10~ 2 , but was hardly worth the effort on vector computers (Cray 
IS), because it makes vectorization difficult and results in short loops. One probably needs 
thousands of vortices for the lumping to bring a significant improvement (the gather-scatter 
capabilities of the newest machines may make a difference). Greengard and Rokhlin report 
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an impressive advantage starting with a few hundred vortices and with e % 10 -5 , again 
on a serial machine (a VAX 8600). The performance of their method on a vector machine 
will be extremely interesting. The nested method is very promising, and the fact the grid 
does not interfere with the body is a definite advantage over vortex-in-cell methods. Also, 
empty cells do not consume any work. The extension to 3-D will require more work since 
we are not just considering the function 1 jz, but is possible. 
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